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1.  Introduction 


Non-crystalline  (amorphous)  ceramics  or  ceramic  glasses  are  used  in  a  variety  of  vital  Army 
personnel,  ground,  and  air  vehicle  applications  that  require  transparent  armor — it  is  ubiquitous  in 
tactical  vehicular  windshields  and  side  windows.  Ceramic  glass  is  inexpensive  and  formable  into 
large,  flat  plate,  and  curved  shapes.  For  many  years  it  has  been  known  that  the  properties  of  glass 
can  be  modified  and  enhanced  through  compositional  modification,  chemical  strengthening, 
annealing,  and  process  control  of  melt  cooling.  Glass  ceramics,  the  controlled  crystallization  of 
nanosized  single  crystals  in  a  glass  matrix,  offer  another  avenue  for  designed  and  enhanced 
property  modifications  for  transparent  and  opaque  material  applications.  In  addition,  certain 
glass  formulations  have  been  shown  to  exhibit  enhanced  performance  against  shaped-charge  jets 
(SCJs)  (figure  1)  and  other  ballistic  threats,  but  it  is  not  understood  why.  This  is  in  part  due  to 
various  short-  and  longer-range  atomic  structural  characteristics  including  atomic  free  volumes, 
cation  coordinations,  bridging  and  non-bridging  oxygen  (O)  atoms,  bonding  energies,  and 
nanoscale  order  characteristics  (short  and  longer  range)  that  are  difficult  or  impossible  to  quantify 
experimentally  for  ceramic  glass. 

In  contrast,  crystalline  ceramics  like  silicon  carbide  (. SiC ),  aluminum  oxynitride  ( AlON ),  and 
others  have  easily  characterizable  micro  structure  s/meso  structures,  which  consist  of  assemblies  of 
individual  single-crystal  grains.  Ceramic  glasses,  on  the  other  hand,  do  not  have  a  conventional 
micro-  or  mesostructure,  as  it  is  understood  for  crystalline  ceramics.  However,  there  are 
microstructural-scale  variations  in  ceramic  glass  that  may  include  density  variations  from  atomic 
free  volume  variations  or  microporosity,  size  of  local  atomic  order,  defects  (inclusions,  large 
pores,  etc.),  and  others  yet  to  be  determined.  The  interaction  of  a  stress/shock  wave  from  a 
dynamic  impact  involves  many  structural  changes  not  easily  characterized  by  conventional 
equations  of  state  and  can  involve  reversible  and  irreversible  densification  and  changes  in  bulk 
short-range  order  structures  comparable  to  phase  changes  in  crystalline  ceramics.  For  example, 
in  simple  Hertzian  indentation  testing,  a  wide  range  of  plastic  or  inelastic  deformation 
mechanisms  have  been  observed  in  a  variety  of  glasses.  Multiscale  computational  design 
methodologies  (figure  2)  for  this  class  of  materials  will,  nevertheless,  require  quantitative  and 
possibly  statistically  based  descriptions  of  the  mesoscale,  although  current  efforts  to  develop  such 
models  have  fallen  far  short  of  this  goal. 
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(a) 


(b) 


Figure  1.  Enhanced  performance  of  SCJs  into  glass  (a)  test  configuration  for  glass  targets, 
and  (b)  penetration  vs.  time  for  several  targets,  after  Moran  et  al.  (7). 


MS  failure  model 


•  Characterization  of  spatial  atomic/microstructural  scale 
glass  structure  and  defect  characteristics 

•  Capturing  these  data  for  embedding  in  meso/macro  code 


(0.1  x  10 


(6  x  10-9m) 


Figure  2.  A  multiscale  model  for  non-crystalline  ceramics  (glass). 
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Our  specific  long-term  research  goals  are  threefold: 


1.  Develop  molecular  dynamics  (MD)  process  models  for  a  series  of  chemically  substituted 
amorphous  silica  (a -Si02  or  fused  silica)  materials  for  the  prediction  of  glass  elastic 
properties  assuming  completely  uniform  glass  “mesostructures.”  If  successful,  such  models 
will  enable  ab  initio  prediction  of  structure-property  relations  in  glass  that  will  be  validated 
with  experimentally  determined  elastic  properties. 

2.  Extend  the  MD  models  to  study  densification  of  chemically  substituted  a -Si02  materials 
under  high  pressures  (~60  GPa  Materials  in  Extreme  Dynamic  Environments  [MEDE]) 
relevant  to  ballistic  events  where  reversible  and  irreversible  density  changes  and  structural 
transformations  have  been  observed.  If  successful,  such  models  will  enable  ab  initio 
prediction  of  a -Si02  compressibility,  kinetics,  and  “glass”  phase  transformations  that  will  be 
used  to  develop  equations  of  state  for  a-Si02  materials,  and  thus  form  a  direct  link  to  the 
continuum  scale. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code 
incorporating  the  effects  of  reversible  and  irreversible  densification,  inelastic  deformation,  and 
overlain  by  a  spatiotemporally  evolving  population  of  growing  defects,  which  coalesce  and 
ultimately  lead  to  fracture  and  fragmentation.  It  is  envisioned  that  at  some  time  in  the  not  too 
distant  future,  fully  concurrent  multiscale  computational  finite  element  codes  will  be  used  by 
analysts  on  a  regular  basis  for  optimal  material  design. 

1.1  Organization  of  the  Report 

The  remainder  of  the  report  is  organized  as  follows.  General  program  objectives  are  outlined  in 
section  2,  and  the  approach  for  modeling  the  multiscale  behavior  of  glass  appears  in  section  3. 
Experimental  work  on  various  glasses  is  described  in  section  4,  which  is  highlighted  by 
indentation  experiments,  edge-on-impact  ballistic  experiments,  ballistic  impact  and  fragmentation 
studies  conducted  at  the  Ernst-Mach  Institute,  and  high  pressure  diamond  anvil  cell  experiments 
conducted  at  the  U.S.  Army  Research  Laboratory  (ARL).  Section  5  describes  the  development  of 
a  new  short-range  pairwise  potential  for  silica  using  ab  initio  molecular  dynamics  (AIMD) 
methods;  the  discovery  of  an  O  —  O  soft  repulsive  shoulder  in  silica  may  explain  multi- stability 
behavior  of  glasses  under  pressure  and  anomalous  densification  behavior.  The  new  potential  is 
used  to  simulate  glass  nanoindentation  experiments,  described  earlier  in  section  4.2. 1 .  Section  6 
introduces  similarity  analysis  for  elastic  media  that  exhibit  fracture  induced  by  indentation;  a 
universal  scaling  law  is  derived,  which  is  invaluable  for  validating  computational  methods 


3 


capable  of  simulating  crack  propagation.  The  mechanics  of  indentation  is  studied  in  section  7, 
where  fully  three-dimensional  stress  and  displacement  fields  are  presented,  which  are  also  useful 
for  verification  of  computational  simulation  methods.  First  principles  quantum  mechanical 
methods  are  used  to  model  densification  and  bulk  modulus  variation  with  pressure  in  section  8. 
Since  inelastic  deformation  and  a  polyamorphic  transformation  occurs  simultaneously  in  fused 
silica,  a  model  is  developed  to  account  for  this  coupled  behavior  in  section  9  and  compared  to 
plate  impact  experiments.  A  peridynamics  computational  code  (section  10)  has  been  developed 
to  enable  modeling  of  quasistatic  and  dynamic  fracture  observed  in  glasses.  Initial  efforts  to 
model  indentation  experiments  and  validation  with  the  universal  scaling  law  are  described  in 
section  6.  Section  11  outlines  a  one-day  short  course  on  the  “Fundamentals  of  Glass  Science,” 
that  was  taught  by  Professor  Arun  K.  Varshneya,  Alfred  University,  at  ARL  on  October  29,  2010. 
Course  attendees  received  copies  of  Varshneya’s  Fundamentals  of  Inorganic  Glasses  (26).  A 
short-term  conceptual  project  to  determine  an  effective  experimental  and  theoretical  approach  to 
model  and  characterize  the  role  of  glassy  materials  in  resisting  ballistic  impact  was  conducted  by 
Professor  Richard  Lehman,  Rutgers  University,  detailed  in  section  12.  Section  13  describes 
ARL’s  new  glass  processing  facility  and  our  initial  efforts  in  processing  borosilicate  systems  for 
evaluation  in  the  transition  of  this  program  to  the  Weapons  and  Materials  Research  Directorate’s 
(WMRD)  core  mission  in  fiscal  year  2013  and  beyond.  Section  14  summarizes  the  conclusions 
of  this  final  report.  Metrics  including  presentations,  publications,  hires,  and  transition  of  the  DSI 
program  to  a  core  mission  program  within  WMRD  are  listed  in  section  15. 


2.  Program  Objectives 


The  long-term  research  goal  of  the  program  is  to  develop  a  concurrent  multiscale  computational 
finite  element  code  for  optimizing  or  enhancing  the  performance  of  various  glasses  against  SCJs; 
the  initial  work  focuses  on  pure  a-Si02  and  chemically  varied  a-Si02  materials.  As  such,  this 
objective  falls  squarely  within  the  purview  of  the  WMRD,  since  multiscale  models  are 
constitutive  models  (specific  to  a  particular  material)  wherein  time-evolving  microstructural 
changes,  such  as  microcrack  growth,  are  fully  coupled  to  the  macroscale,  a  phenomenon  that 
cannot  be  modeled  or  accounted  for  using  classical  homogenization  methods.  A  more  immediate 
research  objective  is  to  understand  why  certain  chemically  substituted  a -Si02  materials  exhibit 
enhanced  performance  in  the  defeat  of  SCJ  and  other  ballistic  threats. 

Our  program  objectives  are  threefold: 
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1.  Develop  MD  process  models  for  a  series  of  chemically  substituted  a-Si02  materials  to  predict 
glass  elastic  properties.  This  glass  plays  an  important  role  in  many  technological  applications 
and  its  structure  has  been  inferred  from  neutron-diffraction,  nuclear  magnetic  resonance,  and 
small-angle  x-ray  scattering  (SAXS)  analysis  to  reveal  a  three-dimensional  network 
consisting  of  tetrahedrally  coordinated  silicon  (Si)  whose  structure  is  constant  throughout  the 
glass  and  defines  its  short-range  order  (SRO).  Long-range  disorder  in  the  structure  is 
manifested  by  a  seemingly  random  variation  in  the  Si-O-Si  bond  angle  in  adjacent  tetrahedra. 
Despite  the  intense  study  of  a -Si02  glass  over  the  last  several  decades,  much  controversy  still 
exists  on  the  best  method  to  model  (i.e.,  via  density  functional  theory,  MD,  Monte  Carlo 
methods,  or  master  equation  techniques)  this  archetypal  material  to  predict  of  elastic 
properties,  diffusivity,  surface  interactions,  bond  angle  distribution,  polyamorphism,  and  melt 
solidification.  Current  models  in  the  literature  are  often  not  fully  validated  and  progress 
towards  this  goal  will  be  made  when  model  predictions  of  elastic  constants  for  a  series  of 
chemically  substituted  a-Si02  glasses  agree  with  experimentally  determined  constants. 

2.  Extend  the  MD  models  to  study  densification  of  chemically  substituted  a-Si02  materials 
under  high  pressures.  Since  long-range  order  in  glass  is  non-existent,  variations  in  the  SRO 
and  intermediate  range  order  (IRO)  must  be  responsible  for  the  enhanced  performance 
observed  in  ballistic  tests  on  certain  a -Si02  glasses.  If  this  is  the  case,  it  may  be  possible  to 
use  MD  models  to  predict  macroscopic  ballistic  performance.  Since  glass  is  subjected  to 
extreme  pressure  and  temperature  during  an  SCJ  event,  it  will  be  necessary  to  study  the 
relationship  between  compressibility,  kinetics,  and  phase  transitions  during  high-pressure 
densification  of  a -Si02  glasses  as  manifested  by  changes  in  coordination  number,  ring  size, 
and  free  volume.  Progress  towards  this  goal  will  be  made  when  MD-derived  equations  of 
state  (EOSs)  agree  with  those  obtained  experimentally  via  diamond  anvil  press  and  plate 
impact  experiments. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code  that 
incorporates  the  effects  of  reversible  and  irreversible  densification  and  inelastic  deformation, 
overlain  with  a  spatiotemporally  evolving  population  of  defects  that  grow,  coalesce,  and 
ultimately  lead  to  fragmentation.  This  objective  will  develop  a  computational  framework  to 
combine  the  objectives  from  (1)  and  (2),  and  incorporate  the  influence  of  fracture  initiation, 
growth,  coalescence,  and  fragmentation  of  surface  and  volume  defects  in  glass  into  a 
comprehensive  concurrent  multiscale  finite  element  model  and  code.  Microcrack  initiation, 
growth,  and  coalescence  (sometimes  referred  to  as  failure  waves)  is  a  multiscale  phenomenon 
that  bridges  all  scales  in  a -Si02  glasses  despite  the  apparent  absence  of  a  structural  mesoscale 
for  this  class  of  materials  (figure  2).  Algorithms  to  develop  fully  two-way  coupled  multiscale 
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codes  are  in  their  infancy,  and  progress  on  this  objective  will  be  realized  with  the  successful 
development  and  implementation  of  a  consistent  scheme  for  coarse-graining  localization 
phenomena  such  as  fracture  failure  observed  in  glass. 


3.  Planned  Approach 


The  planned  approach  consists  of  three  components,  which  are  outlined  in  figure  3: 


1.  Validate  the  MD  models  for  a  series  of  chemically  substituted  a-Si02  materials  to  predict 
glass  elastic  properties.  Although  there  is  no  effort  within  WMRD  to  predict  a -Si02  elastic 
properties,  a  hierarchical  multiscale  modeling  effort  is  currently  underway,  which  is  focused 
on  the  study  of  polycrystalline  (~200  pm  grain  size)  AlON  and  validation  of  quantum  and 
MD  predictions  of  anisotropic  elastic  constants  using  diamond  anvil  cell  (DAC)  and 
focused-ion-beam  (FIB)/scanning  electron  microscopy  (SEM)  compression  tests  on  oriented 
AlON  single  crystals  (27).  We  plan  to  use  MD  methods  (with  possible  MD  coarse-graining) 
to  simulate  glass  process  modeling  during  melt  solidification  by  quenching  a  high-density, 
high  temperature,  and  pressure  (with  possible  polyamorphic  phases)  melt  for  a  series  of 
chemically  substituted  a -Si02  materials.  Next,  the  resulting  room-temperature,  chemically 
modified  structures  will  be  reversibly  deformed  to  predict  the  elastic  properties  that  will  be 
validated  with  experimentally  determined  elastic  properties. 

2.  Validate  the  MD  models  for  densification  of  chemically  substituted  a-Si02  materials  under 
high  pressures.  Although  there  is  currently  no  effort  within  WMRD  to  predict  the  EOS  of 
chemically  substituted  a -Si02  materials,  MD  methods  have  been  used  to  predict  high- 
pressure  densification  in  these  materials.  MD  simulations  of  pure  a -Si02  materials  reveal  a 
Hugoniot  elastic  limit  (HEL)  of  about  10  GPa,  and  an  anomalous  maximum  in 
compressibility  at  around  3  GPa.  Experiments  where  samples  have  been  compressed  to 
pressures  lower  than  10  GPa  are  indistinguishable  from  the  original  material,  whereas  above 
10  GPa,  materials  can  sustain  an  irreversible  density  increase  from  10%-20%  higher  than  the 
starting  material,  although  there  is  controversy  as  to  whether  the  mechanism  is  due  to 
irreversible  coordination  defects  or  permanent  ring  size  modification.  In  contrast  to  the 
behavior  of  a -Si02,  crystalline  quartz  (a-Si02)  undergoes  very  well-known  high  pressure 
polymorphic  phase  transitions  into  a  Coesite  phase  and  a  Stishovite  phase  (figure  4),  which 
involve  changes  in  coordination  of  the  Si  cation  from  four  to  six  O  atoms.  A  combination  of 
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Figure  3.  Five-year  roadmap  consistent  with  the  WMRD  brittle  materials  program. 


diamond  anvil  press  and  plate  impact  experiments  will  be  conducted  on  a  series  of  chemically 
substituted  a-S702  materials  and  compared  with  MD  densification  simulations  in  glass  in 
order  to  understand  the  influence  of  glass  modifiers  on  changes  in  the  shock  response  of  these 
materials.  EOSs  for  a  subset  of  promising  chemically  substituted  materials  will  be  developed 
and  implemented  into  a  continuum  code  to  determine  if  any  of  the  chemically  substituted 
materials  exhibit  enhanced  ballistic  performance. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code  that 

incorporates  the  effects  of  reversible  and  irreversible  densification  and  inelastic  deformation 
overlain  with  a  spatiotemporally  varying  population  of  defects  that  grow,  coalesce,  and 
ultimately  lead  to  fragmentation.  The  ultimate  research  objective  is  to  develop  a 
physics-based  multiscale  computational  finite  element  code  for  studying  densification  and 
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T  (°C) 


Figure  4.  Crystal  structures  of  quartz,  after  Frye  (2). 


dynamic  fracture  in  non-crystalline  ceramics  (see  figure  2).  Atomistic  behavior  will  be  linked 
to  macroscopic  elastic  properties  and  densification  behavior  through  development  of  an  EOS 
from  first  principles  as  outlined  in  components  (1)  and  (2).  At  this  stage,  what  remains  to  be 
accomplished,  is  to  successfully  link,  in  a  concurrent  fashion,  multiscale  failure  phenomena  in 
a-Si02  materials  by  incorporating  the  important  role  that  pre-existing  surface  and  volume 
defects  have  on  the  microcrack  growth,  coalescence,  and  fragmentation  in  this  class  of 
materials.  Over  the  past  five  years,  the  first  author  has  also  been  directly  involved  in 
development  of  a  parallel,  concurrent  multiscale  code  for  heterogeneous  viscoelastic 
composites  (28)  under  the  auspices  of  an  ARL/University  of  Nebraska  cooperative 
agreement,  which  will  be  leveraged  and  used  as  the  framework  for  the  development  of  a 
concurrent  multiscale  model  of  a-Si02  materials. 


The  chief  challenge  for  brittle  materials  is  to  correctly  account  for  the  growth  kinetics  of 
microcracks  in  a  multiscale  computational  environment.  The  propagation  of  free  internal 
boundaries  at  lower  scales  will  be  “coarse-grained”  to  higher  scales,  where  global  fracture 
failure  and  fragmentation  is  observed.  As  such,  coarse-graining  algorithms  will  need  to  be 
validated  through  continuum-scale  experiments  on  a -Si02  materials  that  measure  dynamic 
crack  propagation  speeds,  mixed-mode  failure,  and  crack  bifurcation  phenomena  using 
coherent  gradient  sensing  and  high-speed  imaging  techniques;  ARL  possesses  capabilities  for 
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conducting  such  dynamic  fracture  experiments  in  a-SiO-2  materials  through  the  recent 
establishment  of  a  coherent  gradient  sensing/imaging  facility  funded  by  the  ongoing 
multiscale  modeling  effort  of  AlON .  Models  and  validation  of  the  initiation  and  propagation 
of  discrete  fractures  in  a-SiC>2  materials  should  transition  naturally  into  models  of 
fragmentation  and  comminution  for  behind-armor-debris  applications.  Fragmentation 
experiments  have  classically  been  conducted  using  dynamically  expanding  ring  experiments 
for  defining  fragment  size  versus  strain  rate  and  will  be  used  to  validate  computational  models 
of  fragmentation.  The  development  of  consistent  coarse-graining  algorithms  for  fracture  in 
materials,  which  is  associated  with  failure  and  loss  of  material  stability,  is  largely  unexplored 
and  is  the  primary  high-risk  goal  of  this  section. 


4.  Experimental  Work  and  Background 


4.1  Background 

4.1.1  Compositions 

Silicate-based  ceramic  glasses  are  based  on  chemical  substitutions  into  a  SiO 4  tetrahedral-based 
polymeric-like  structure;  fused  silica  is  an  amorphous  (non-crystalline)  form  of  pure  SiO 2.  Table 
1  lists  the  compositions  and  properties  of  typical  glasses. 


Table  1 .  Glass  compositions  in  %  and  selected  properties. 


Glass 

Si02 

Aho^ 

CaO 

B2O3 

Ncl20 

k2o 

p  (g/cm3) 

E  (GPa) 

V 

Borofloat 

80.5 

2.5 

0.02 

12.7 

3.55 

.64 

2.23 

62.3 

.207 

Starphire 

73.2 

1.44 

10.27 

- 

14.72 

.01 

2.51 

72.3 

.23 

Fused  Silica 

100 

- 

- 

- 

- 

- 

2.2 

73.0 

.17 

4.1.2  Structural  Characteristics  of  Ceramic  Glasses 

Simplistically,  the  predominant  macro-characteristics  (micron  and  larger)  can  be  a  variety  of 
defects  including  inclusions,  bubbles,  large  pores,  and  residual  compressive  or  tensile  stress.  The 
notion  of  an  array  of  crystalline  grains  separated  by  grain  boundaries  (a  microstructure)  does  not 
exist  in  glass.  Rather,  there  is  a  complete  lack  of  long-range  order,  but  short-  and 
intermediate-range  order  at  the  nano  structural  scale: 
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Short-range  order:  Mostly  atom  to  atom  bond  lengths,  less  than  0.5  nm  and  bond  angles; 
characterized  by  radial  distribution  functions  (RDFs)  and  SAXS. 


•  Intermediate-range  order:  In  silica-based  glasses,  this  is  the  polymerization  of  the  silica 
atomic  tetrahedra  (one  Si  atom  surrounded  by  four  oxygen  atoms)  into  various  size  ring 
structures  of  joined  tetrahedra,  which  can  consist  of  4,  5,  6,  7,  8,  or  so  ring  groups  of 
tetrahedra.  Substitution  of  other  cations  (Na,  K,  Mg,  Ca,  etc.)  and  B  into  silica-based 
glasses  can  have  profound  effects  on  the  IRO.  It  is  also  important  to  note  that  B  bonds  to  three 
oxygen  atoms  in  flat  triangles. 

•  Free-volume:  In  crystalline  materials,  using  the  theoretical  density,  it  is  straightforward  to 
calculate  the  atomically  unoccupied  free  space.  In  glasses,  this  unoccupied  atomic  space  is 
referred  to  as  “free  volume,”  but  because  of  an  unknown  periodic  structure,  it  is  extremely 
difficult  to  quantify  in  glasses.  The  free  volume  plays  a  critical  role  in  glass  densification 
under  stress/pressure. 

•  There  can  be  significant  complex  spatial  variations  of  defects,  free  volume,  atomic  structure, 
and  resulting  properties  at  the  nanoscale. 

•  It  is  the  current  wisdom  of  the  glass  community  that  the  short-  and  intermediate -range  order  at 
the  nanoscale  in  glasses  can  have  a  significant  influence  on  some  properties. 

4.1.3  Effect  of  Stress/Pressure 

Deformation  and  failure  in  ceramic  glasses  begins  with  reversible  to  irreversible  densification 
and/or  critical  cracks  nucleating  at  defects: 


•  High  pressure  can  have  significant  effect  on  coordination  of  Si  changing  from  typical  fourfold 
to  sixfold  coordination  of  oxygen  around  Si  -  (quartz-like  to  coesite-like  to  stishovite-like 
structures,  as  for  crystalline  quartz  shown  in  figure  4). 

•  The  degree  of  polymerization  or  distribution  of  the  ring  structures  can  also  change  as  a 
function  of  pressure. 

•  The  common  consensus  is  that  the  ring  distribution  is  the  primary  control  of  some  properties; 
however,  at  very  high  pressures,  the  change  from  four  to  sixfold  Si  coordination  will 
significantly  influence  properties  as  well. 


10 


•  Some  properties  of  fused  silica  (pure  SiO 2  glass)  are  anomalous,  e.g.,  a  negative  change  in 
shear  modulus  as  a  function  of  pressure.  Bando  et  al.  (3)  show  that  the  radius  of  curvature  of 
a  crack  in  glass  (figure  5)  can  be  about  1.5  nm,  suggesting  that  the  IRO  can  significantly 
influence  crack  propagation. 


Figure  5.  The  radius  of  curvature  of  a  crack  in  glass  after  Bando  (3). 
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4.1.4  Plasticity  in  Glass 


Ito  ( 4 )  has  presented  fairly  simple  methods  of  determining  the  brittleness  (figure  6)  or,  conversely, 
the  plasticity  of  glasses,  which  he  uses  to  suggest  that  the  brittleness  is  dependent  of  IRO  or  the 
distribution  of  ring  structures  seen  in  figure  7.  Table  2  lists  values  for  the  brittleness  parameters. 


Figure  6.  Brittleness  vs.  density  for  glasses  in  the  S1O2  and  -based 
glasses,  after  Ito  ( 4 ). 


(a)  (b) 

Figure  7.  (a)  Structure  of  soda  lime  glass  by  MD  simulation  where  the 

number  shown  is  the  ring  size  and  (b)  Number  of  network  rings 
in  soda  lime  (SL)  and  less  brittle  (LB)  glasses;  LB  glass  are 
more  polymerized  than  SL  glass  after  Ito  (4). 
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Table  2.  Brittleness  parameters. 


Glass 

Brittleness  Parameter  (gm  l/2) 

Fused  Silica 

~10 

B2Os  Based  glass 

rsJ  1 

Soda  Lime 

~  5-7 

The  brittleness,  therefore,  seems  to  be  dependent  on  the  deformation  and  fracture  behavior,  which 
depends  on  flow  and  densification  before  crack  initiation  and  on  the  bond  strength  of  the  network 
and  seems  to  decrease  with  a  decrease  in  density — a  free  volume  effect.  This  is  addressed  by 
Ito  (4)  in  the  same  paper. 

Conclusions  from  this  work  are  as  follows: 


•  Both  glasses  are  commercial  SL  based  glasses:  (. Na ,  K)20  —  (Mg,  Ca)0  -  Si02. 

•  The  LB  glass  appears  to  have  a  higher  polymerized  network  than  the  SL  glass,  i.e.,  a 
significant  difference  in  the  ring  structure  distribution. 

•  IRO  at  the  nanoscale  seems  to  be  controlling  the  brittleness  of  these  glasses. 

4.2  Experimental  Results 

The  absorption/dissipation  of  deposited  energy  in  an  extreme  impact  event  depends  on  the  various 
deformation  and  failure/fracture  mechanisms  that  are  activated  during  the  event.  In  addition,  in  a 
multiscale  modeling  and  simulation  “Protection  Materials  by  Design”  approach,  it  is  absolutely 
critical  to  experimentally  determine  the  key  properties  at  the  various  scales  to  validate  the 
theoretical  computational  results.  We  have  used  various  quasi-static  indentation,  edge-on-impact 
(EOI),  and  ballistic  impact  tests  for  this  purpose. 

4.2.1  Indentation 

Studies  on  the  deformation  and  fracture  of  glasses  and  AlON  using  a  spherical  500-yum-diameter 
diamond  indenter  was  recently  carried  out  by  Wilantewicz  (5).  Results  for  a  SL  (Starphire), 
boron  substituted  glass  (Borofloat),  fused  silica,  and  AlON  are  illustrated  in  figure  8  and  tables  3 
and  4.  There  are  significant  differences  in  the  deformation  and  fracture  behavior  of  these 
materials. 
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2  5  10  20  30  35  40  45  50  65  6  5  75  100  125  150  200 

Load  (N) 


Elastic  recovery  and  percent  elastic  recovery  as  a  function  of  the  maximum  indentation  load  for  AiON,  Starphire 
fused  silica  and  Borofioat glasses  using  a  500  mm  diameter  spherical  diamond  indenter,  ten  indentations  at 
each  load  were  made.  "Failure  Behavior  of  Glass  and  Aluminum  Oxynitride  {AION) 

Tiles  Under  Spherical  Ind  enters”,  Trevor  E.  Wilantewicz.  ARL-TR45180  May  2010 


Figure  8.  Elastic  recovery  vs.  load  for  a  variety  of  glasses  after  Wilantewicz  (5). 


Table  3.  Elastic  recovery  in 

indentation  tests  at  200  N. 


Material 

%  Elastic  Recovery 

AION 

71 

Starphire 

73 

Borofioat 

79 

Fused  silica 

86 

Table  4.  Deformation  and  fracture  loads  after  Wilantewicz  (5). 


Onset 

All  Tests 

Onset  Ring 

All  Tests 

Onset  Radial 

All  Tests 

Dimpling 

Dimpled 

Cracking 

Ring  Cracked 

Cracking 

Radial  Cracked 

Material 

(N) 

(N) 

(N) 

(N) 

(N) 

(N) 

Starphire  (tin) 

30 

30 

65 

75 

75 

100 

Starphire  (air-) 

20 

30 

65 

100 

100 

125 

Borofioat  (tin) 

30 

35 

30 

45 

100 

200 

Silica  Glass 

75 

100 

20 

30 

65 

75 

AION 

35 

45 

45 

65 

40 

75 
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It  is  clear  that  these  materials  behave  in  significantly  different  ways.  The  onset  of  dimpling  is  the 
result  of  a  permanent  plastic  deformation.  The  normal  expectation  for  these  materials  is  that  as  a 
function  of  increasing  indentation  load  the  material  response  would  proceed  through  an  elastic 
regime,  then  through  a  plastic  regime,  and  finally  into  a  cracking/fracturing  regime.  Table  5 
summarizes  the  observations.  Silica  glass  (fused  silica),  however,  behaves  in  a  dramatically 
different  way,  fracturing  prior  to  a  plastic  mechanism. 


Table  5.  Onset  of  elastic,  plastic,  and 
fracture  responses  after 
Wilantewicz  (5). 


Material 

Response 

Starphire: 

elastic  — >•  plastic  — >•  fracture 

Borofloat: 

elastic  — >•  plastic  — >•  fracture 

Silica  Glass: 

elastic  — >•  fracture  — >•  plastic 

AlON: 

elastic  — >•  plastic  — »  fracture 

4.2.2  Edge-on-impact  Studies  of  Fused  Silica 

The  experimental  arrangement  is  illustrated  in  figure  9  and  a  series  of  time-resolved  photographs 
are  presented  in  figure  10. 


Plane 


Figure  9.  EOI  experimental  arrangement  after  Strassburger  et  al.  (6). 
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Test  No  14878  t  =  8.7  ps  Fig.  4:  Test  No  14878:  t»  9.7  ps  «=10.7ps 

Measurement  of  transversal  Measurement  of  crack  retoaty  from  crack  center  expansion 

wave  speed  and  crack  velocity 


(a) 


(b) 


Figure  10.  (a)  Illustration  of  a  series  of  EOI  tests  in  fused  silica  by  a  solid  steel 
cylinder  at  350  m/s  at  various  times;  first  and  third  rows  illustrate 
shadowgraph  photographs  in  plane  light  showing  damage;  second  and 
fourth  rows  arc  photos  in  crossed  polarized  light,  which  show  the 
propagation  of  stress  via  a  photoelastic  effect;  and  (b)  illustrates  the 
irregular  nature  of  the  damage  front  due  to  the  presence  of 
macro-defects  from  the  same  test  after  Strassburger  et  al.  (6). 
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Table  6  lists  the  measured  velocities  of  the  longitudinal  waves,  transverse  waves  (shear  waves), 
and  crack  and  damage  front  velocities  in  fused  silica.  Note  that  the  longitudinal  wave  velocity 
for  fused  silica  is  5.93  km/s  and  the  shear  wave  velocity  is  3.77  km/s. 


Table  6.  Compilation  of  measured  wave,  crack,  and  damage  velocities  in 
fused  silica  after  Strassburger  et  al.  (6). 


Impact  velocity  (m/s) 

150 

260 

350 

Longitudinal  wave  speed  (m/s) 

shadowgraphs 

5975 

6076 

5823 

Longitudinal  wave  speed  (m/s) 

crossed  polarizers 

5814 

5796 

5491 

Transverse  wave  speed  (m/s) 

shadowgraphs 

- 

3500 

3670 

Crack  velocity  (m/s) 

shadowgraphs 

2234 

2149 

2120 

Damage  velocity  (m/s) 

shadowgraphs 

5641 

5728 

5121 

4.2.3  Visualization  and  Analysis  of  Ballistic  Impact  Damage  and  Fragmentation  in 
Various  Glass  Plates 

In  this  series  of  experiments  Borofioat,  Starphire,  and  fused  silica  were  tested  in  a  standard 
ballistic  configuration.  The  plates  were  impacted  by  a  7.62-mm  armor-piercing  (AP)  round,  and 
a  solid  steel  cylinder  inside  of  a  box  to  contain  all  of  the  resulting  fragments.  The  fragments 
were  removed  from  the  box  with  a  vacuum  and  then  sorted  by  sieves.  Figure  1 1  illustrates  the 
experimental  arrangement. 

Table  7  lists  the  details  of  the  various  tests  conducted  on  the  three  glass  types.  Note  that  the 
dimensions  of  the  Borofioat  glass  plate  used  in  test,  #17742,  were  significantly  different  than  the 
others,  which  has  skewed  the  fragmentation  results  at  the  largest  sieve  size  2  mm.  Figure  12 
illustrates  a  series  of  very-high-speed  photographs  as  a  function  of  time  for  the  four  tests  listed  in 
table  7. 


Table  7.  Test  parameters  with  different  types  of  glass. 


EMI  Test  No. 

Type 

Dimensions  (mm) 

Thickness  (mm) 

Projectile 

Impact  Velocity  (m/s) 

17742 

Borofioat 

149.4  x  149.7 

12.94 

cylinder 

1089 

17749 

Fused  silica 

101.65  x  101.67 

12.75 

cylinder 

1107 

17750 

Fused  silica 

101.62  x  101.62 

12.75 

7.62  AP 

824 

17751 

Starphire 

99.9  x  99.7 

10.06 

cylinder 

1115 

17 


(a) 


(b) 


Figure  1 1 .  Schematic  of  (a)  ballistic  test  configuration  and  (b)  target  after 
Strassburger  et  al.  (7). 


#  17751:  soda  lime:  vp  =  #  17742:  Borofloat. 

1115m/s:Solid  cylinder  solid  cylinder:  1089m/s 


#  17749:fused  silica: 
solid  cylinder:  1100  m/s 


■HU  11 


K  X  15  vs 


#  17750:fused  silica: 
7.62  AP:  800  m/s 


Figure  12.  Selection  of  high-speed  photographs  from  impact  on  various  glasses 
after  Strassburger  et  al.  (7). 
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The  observed  propagation  velocities  of  the  damage  zone  under  impact  of  the  steel  cylinder  are  all 
below  the  transverse  wave  velocities  of  the  materials,  as  seen  in  table  8,  which  compiles  the  wave 
and  fracture  velocity  data  determined  from  EOI  tests  along  with  data  from  the  actual  test  series. 
When  fused  silica  was  impacted  by  the  7.62-mm  AP  projectile,  the  formation  of  single  radial 
cracks  were  observed,  which  propagated  at  an  average  velocity  of  2394  m/s,  which  is  in  very 
good  agreement  with  the  crack  velocity  determined  for  fused  silica  from  EOI  tests. 


Table  8.  Compilation  of  wave  and  fracture  velocity  data  after  Strassburger  et  al.  (7). 


Glass  Type 

Longitudinal  Wave 
Velocity  (m/s) 

Transverse  Wave 

Velocity  (m/s) 

Terminal  Crack 

Velocity  (m/s) 

Damage 
Velocity  (m/s) 

Starphire 

5890 

3570 

1580 

3073® 

Borosilicate 

5543 

- 

2034 

2857® 

Fused  Silica 

6021 

3858 

2400 

3007®  23941 

®  steel  cylinder,  v  =  1100  m/s 
1  AP  projectile,  v  =  824  m/s 


4.2.4  High-speed  Photography 

The  three  types  of  glass  tested  exhibited  a  very  similar  fracture  pattern  under  impact  of  a  steel 
cylinder  at  1100  m/s.  A  circular  damage  zone  developed  around  the  impact  site,  in  which  the 
glass  was  strongly  fragmented  and  no  single  cracks  could  be  distinguished  during  the  first  10  fis 
after  impact.  After  this  first  damage  propagation  phase,  single  cracks  became  discernible  at  the 
perimeter  of  the  damage  zone.  In  the  following,  a  severely  damaged  inner  zone  and  an  outer 
zone  with  a  lower  fracture  density  could  be  distinguished.  Selections  of  16  high-speed 
photographs  from  each  test  with  the  cylindrical  projectile  are  presented  in  figures  13,  15,  and  17. 
The  analysis  of  damage  propagation  revealed  different  velocities  with  the  three  types  of  glass. 

The  corresponding  path-time  histories  are  depicted  in  figures  14,  16,  and  18. 

When  fused  silica  was  impacted  by  the  7.62-mm  AP  projectile  the  formation  of  single  radial 
cracks  could  be  observed,  which  propagated  at  an  average  velocity  of  2394  m/s.  This  is  in  very 
good  agreement  with  the  crack  velocity  determined  for  fused  silica  from  EOI  tests  (29).  A 
selection  of  16  high-speed  photographs  from  the  AP  projectile  impact  on  fused  silica  at  824  m/s  is 
presented  in  figure  19.  The  corresponding  path-time  history  of  fracture  propagation  shows  the 
diagram  in  figure  20. 
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Figure  13.  Selection  of  high-speed  photographs  from  impact  on  Borofloat  glass  at 
1089  m/s;  Test  #17742  . 


Figure  14.  Position-time  history  of  damage  propagation  in  Borofloat  glass  at  1089 
m/s;  Test  #17742  . 
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Figure  15.  Selection  of  high-speed  photographs  from  impact  on  fused  silica  glass 
at  1 107  m/s;  Test  #17749  . 


Time  [ps] 


Figure  16.  Position-time  history  of  damage  propagation  in  fused  silica  glass  at 
1107  m/s;  Test  #17749  . 
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Figure  17.  Selection  of  high-speed  photographs  from  impact  on  SL  glass  at  1115 
m/s;  Test  #17751  . 


0  5  10  15 

Time  [ps] 


Figure  18.  Position-time  history  of  damage  propagation  in  SL  glass  at  1115  m/s; 
Test  #17751  . 
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Figure  19.  Selection  of  high-speed  photographs  from  impact  on  fused  silica  glass 
at  824  m/s;  Test  #17750  . 


0  5  10  15  20 


Time  [ps] 


Figure  20.  Position-time  history  of  damage  propagation  in  fused  silica  glass  at 
824  m/s;  Test  #17750  . 
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4.2.5  Fragmentation  Analysis  with  Glass 


The  glass  fragments  from  these  experiments  were  collected  and  separated  into  size  classes  by  a 
chain  of  sieves  in  the  same  way  as  the  ceramic  fragments.  The  mesh  sizes  used  were  2  mm,  1 
mm,  0.5  mm,  200  /tm,  100  /im,  63  /mi,  and  25  /fm.  The  total  mass  of  each  size  fraction  was 
determined.  Figure  21  presents  the  values  of  the  total  fragment  mass  in  the  different  size  classes 
with  the  four  different  configurations.  The  corresponding  cumulative  mass  plot  is  shown  in 
figure  22. 

A  very  high  total  mass  of  fragments  of  size  >  2  mm  can  be  recognized  with  the  Borofloat  glass. 
This  result  can  be  attributed  to  the  size  of  the  sample  (150  mm  x  150  mm),  and  therefore,  the 
higher  total  mass  (643  g)  compared  to  the  fused  silica  (290  g)  and  the  SL  glass  samples  (250  g). 
The  highest  fragment  mass  was  found  with  fused  silica  in  the  size  classes  from  63  m  to  1  mm. 
The  cumulative  mass  plot  also  reflects  the  highest  degree  of  fragmentation  with  fused  silica  in  the 
tests  with  the  steel  cylinder.  The  least  degree  of  fragmentation  was  observed  with  fused  silica 
impacted  by  the  AP  projectile.  In  this  case,  the  projectile  also  penetrated  the  aluminum  backing 
nearly  complete  and  a  lower  amount  of  energy  was  dissipated  in  the  interaction  with  the  glass. 


2  1  0.5  0.2  0.1  0.063  0.025 

mesh  size  [mm] 


Figure  21.  Fragment  mass  distribution  from  sieve  analysis  for  tests  with  various 
glasses  after  Strassburger  et  al.  (7). 
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mesh  size  [mm] 


Figure  22.  Cumulative  mass  plot  for  tests  in  various  glasses. 
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4.2.6  Density  Measurements  with  Glass  Fragments 


The  compressibility  of  glass  has  been  studied  for  several  years  and  permanent  densification  up  to 
20%  under  high  pressure  has  been  reported  (8,  30).  Figure  23  shows  the  results  of  densification 
measurements  with  silica  glass  at  25  °C  from  three  different  researchers. 


Figure  23.  Densification  of  silica  glass  at  25  °C  as  a  function  of  pressure  from 
Mackenzie  (8);  curve  A  from  Roy  and  Cohen  (9),  curve  B  from 
Christiansen  et  al.  (10),  and  curve  C  from  Bridgman  (11). 


The  elastic  shock  induced  at  the  impact  surface  between  a  projectile  and  a  target  is  determined  by 
the  acoustic  impedances  of  the  projectile  zp  and  target  material  zp  and  is  given  as  follows  (31): 

P=  ZpZT  Vp  =  PPCrPTCT  yp  (1) 

Zp  +  Zp  PpCp  +  ppCp 

where  pp,  pp  is  the  density,  and  c P  and  cp  are  the  longitudinal  sound  wave  velocity  of  the 
projectile  and  target  material,  respectively. 

If  a  steel  projectile  impacts  a  fused  silica  target,  we  have  pp  =  7.85  g/cm3,  pT  =  2.2  g/cm3,  cp  = 
5100  m/s,  and  cp  =  6021  m/s,  then  a  shock  pressure  of  ~1 1  GPa  (~1 10  Kbar)  is  generated  during 
impact  of  a  steel  projectile  traveling  at  about  Vp  =  1  km/s  onto  a  fused  silica  target.  As  per  the 
data  in  figure  23  (curves  A,  B),  a  4%-10%  densification  could  possibly  be  expected. 
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Therefore,  floatation-sink  tests  were  performed  with  fused  silica  fragments  of  size  1-2  mm  and 
0. 1-0.2  mm  from  impact  tests  #17749  (steel  cylinder,  1107  m/s)  and  #17750  (7.62-mm  AP 
projectile,  824  m/s).  Aqueous  sodium  polytungstate  (NaPW)  was  used  as  heavy  liquid  for  the 
float-sink  analysis.  The  density  of  the  heavy  liquid  was  varied  within  the  range  2.19-2.4  g/ml. 
The  density  was  measured  by  means  of  precision  areometers.  Table  9  summarizes  the  test  data. 


Table  9.  Data  from  floatation-sink  tests  with  fused  silica  fragments  from  impacted  samples. 


Test 

Impact 

Particle 

Mass 

Volume 

Addition 

Addition 

Density 

Observation: 

Test  No. 

Size 

sample 

NaPW-soln: 

dist.  H2O 

NaPW 

NaPW-soln: 

(mm) 

(g) 

(ml) 

(ml) 

(ml) 

(g/ml) 

1. 

HI 

1.000 

4.99 

50 

15 

X 

2.40 

silica  swims 

1. 

1.000 

5.08 

300 

25 

X 

2.35 

silica  swims 

1. 

1.000 

5.26 

325 

15 

X 

2.29 

silica  swims 

1. 

1.000 

5.01 

335 

15 

X 

2.24 

silica  swims 

1. 

HI 

1.000 

5.09 

350 

15 

X 

2.19 

silica  swims 

2. 

17750 

1.000 

5.06 

350 

15 

X 

2.19 

silica  sinks 

3. 

17749 

0.100 

5.35 

350 

15 

X 

2.19 

silica  sinks 

4. 

17750 

0.100 

3.24 

350 

15 

X 

2.19 

silica  sinks 

5. 

17749 

0.100 

2.25 

255 

X 

20 

2.25 

silica  swims 

6. 

17750 

0.100 

2.43 

255 

X 

20 

2.25 

silica  swims 

With  the  analyzed  fragments,  a  densification  could  not  be  measured.  On  one  hand,  the  collected 
fragments  could  not  be  allocated  to  their  original  position  in  the  fused  silica  sample.  Therefore,  it 
was  not  possible  to  determine  whether  the  analyzed  fragments  originated  from  the  impact  zone 
and  had  been  subject  to  high  pressure.  On  the  other  hand,  due  to  the  small  diameter  of  the 
projectile,  the  duration  of  the  pressure  pulse  could  only  have  been  short,  because  of  the  release 
waves  from  the  edge  of  the  projectile. 
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4.2.7  Micro/Nanocrystallinity  in  Glass 


In  discussions  with  Professor  Adrian  Wright  from  Reading  University,  UK,  an  expert  in  the 
structure  of  glass,  he  has  said  that  the  theory  of  glass  structure  is  still  evolving.  The  IRO  is 
important,  but  poorly  understood,  and  the  current  consensus  is  that  there  is  chemical 
nanoheterogeneity,  which  may  lead  to  nanoislands  of  crystallinity.  Thus,  the  nature  of  IRO  in 
glass  is  still  evolving  and  controversial.  Work  by  Saito  et  al.  (32)  using  light  scattering  studies 
on  silica  glass  concluded  that  medium  (intermediate)-range  order  structures,  microcrystallites, 
exist  in  silica  glass.  Using  the  polarization  ratio,  a  measure  of  anisotropy,  suggests  that  the 
depolarized  scattering  is  only  attributed  to  the  density  fluctuations  with  microscopic  anisotropic 
scattering  elements.  Anisotropy  of  such  scattering  elements  could  be  explained  by  the  existence 
of  microcrystallites. 

4.2.8  Thermal  Shock  Behavior  of  Glasses 

Bradt  and  Martens  (33)  in  a  recent  article  have  reviewed  the  thermal  shock  resistance  of  various 
glasses,  including  borosilicate  and  SL  glasses.  They  show  that  thermal  stresses  that  develop 
during  temperature  changes  are  the  primary  cause  of  fracture  initiation  and  failure.  Using  a 
simple  formula  for  the  generation  of  linear  elastic  thermal  stresses  ats. 


ats  =  aEAT 


(2) 


where,  a  is  the  coefficient  of  thermal  expansion,  E  is  Young’s  modulus,  and  AT  is  the  local 
temperature  differential.  Rewriting  equation  2  gives  the  temperature  difference  required  to  cause 
failure  at  a  stress  07, 


AT 


_T/_ 

aE 


(3) 


hence,  using  this  equation  one  can  approximate  the  AT  for  glass  fracture.  Assuming  the 
properties  for  a  typical  SL  glass  and  a  borosilicate  glass  (table  10),  the  AT  for  fracture  can  be 
approximated. 


Other  more  complex  equations  can  be  used,  but  reflect  very  similar  general  trends.  The  general 
conclusion,  therefore,  is  that  the  linear  thermal  expansion  of  glasses  is  a  major  factor  in  initiating 
crack  formation  and  ultimate  failure. 
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Table  10.  AT  required  to  fail  SL  and 
borosilicate  glasses. 


Soda  lime 

Borosilicate 

af 

30  MPa 

30  MPa 

E 

68  GPa 

62  GPa 

a 

9  x  lO-60^”1 

3  x  10"6oC_1 

AT  for  a f 

55  °C 

183°G 

4.2.9  Diamond  Anvil  Cell  (DAC)  Experiments  to  Study  Density  Changes  in  Glass  with 
Pressure 

The  Army  has  a  longstanding  interest  in  understanding  the  properties  of  transparent  protection 
materials,  which  form  the  basis  of  windshields,  optics,  and  viewports.  The  silicate  glass  family 
of  materials  makes  up  the  traditional  choices  for  these  applications.  These  materials  fall  into  a 
much  broader  category  of  materials  called  amorphous  glasses.  The  silicate  glasses  have  several 
subcategories  determined  by  the  secondary  component,  namely,  the  Borofloat  and  SL  varieties 
studied  in  this  project.  The  major  interest  from  the  Army  perspective  in  these  glasses  has  to  do 
with  the  differences  in  ballistic  properties  while  being  compositionally  similar.  The  term 
“amorphous  materials”  when  applied  to  glass  is  slightly  misleading.  These  systems  are  thought 
to  have  several  different  arrangements  from  a  rigid  short-range  silicon-oxygen  bond  length,  to 
intermediate-range  ring  structures,  and  finally,  to  an  overall  amorphous  packing  of  these  smaller 
structures  on  the  mesoscale.  This  structural  hierarchy  is  still  theoretical  as  the  amorphous  nature 
makes  study  and  analysis  difficult  with  existing  techniques.  Furthermore  as  silicate  glasses  have 
identical  distances  for  the  basic  silicon-oxygen  bond,  it  is  apparent  that  the  differences  in  the 
glass  performance  are  due  to  the  intermediate -range  packing  and  ring  structures.  Understanding 
this  difference  is  one  of  the  goals  of  this  research. 

The  current  state  of  the  art  in  understanding  the  IRO  of  silicate  glasses  is  work  involving 
complementary  techniques  of  modeling  and  experiments.  For  the  purposes  of  understanding  how 
these  materials  change  under  ballistic  impact,  static  data  allow  for  comparisons  to  shock 
Hugoniot  data.  There  has  been  previous  work  involving  static  pressure  experiments  on  fused 
silica  (34,  35),  and  also  more  recently  on  the  similar  germanium  glass  (36).  The  techniques 
employed  in  these  studies  were  neutron  diffraction  and  x-ray  diffraction  combined  with  structure 
factor  analysis  on  the  experimental  side  combined  with  MD  simulations  on  the  theoretical  side. 
Even  with  these  detailed  studies,  the  mechanism  for  densification  is  still  elusive,  without  adding 
the  additional  complexity  of  the  effects  of  the  secondary  components. 
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The  DAC  provides  a  convenient  method  of  studying  glass  samples  under  pressure.  The 
experimental  technique  involves  compression  of  the  sample,  the  hydrostatic  pressure  medium, 
and  the  metal  gasket  between  opposed  diamond  anvils.  The  hydrostatic  pressure  media  use  for 
these  experiments  is  argon  gas  and  the  gaskets  were  composed  of  either  rhenium  or  steel.  The 
standard  ruby  fluorescence  pressure  scale  was  used  (37).  Samples  were  loaded  in  the  DAC  and 
then  monitored  for  changes  in  sample  area  with  pressure.  This  area  was  then  converted  to 
volume  measurements  using  the  method  from  previous  work  investigating  the  effect  of  helium  on 
compression  of  Si02  glass  (12).  The  pressure  media  is  assumed  to  be  hydrostatic  in  this 
technique.  The  equation  used  is  reproduced  below: 


V  A  / - 

T7  =  “T  'J  AIA  0  ’  (4) 

ro  A) 

where  V  and  Vo  refer  to  final  and  initial  volumes,  and  A  and  /10  refer  to  final  and  initial  areas. 

The  samples  were  initially  prepared  using  a  polishing  technique  to  achieve  the  required  thickness 
of  <  20  microns.  The  samples  were  then  extracted  yielding  a  glass  chip  that  was  approximately 
30-50  microns  in  diameter.  This  was  so  that  the  sample  could  be  loaded  in  a  gasket  with 
100-micron  hole.  It  was  noted  that  the  irregular  surface  of  the  chips  made  interpretation  of  the 
edges  for  the  area  calculation  difficult.  In  addition,  it  was  thought  that  the  comers  could  provide 
stress  points  in  the  sample  and  introduce  unintended  strain  into  the  measurement.  To  correct  this 
issue,  the  samples  were  machined  using  a  FIB  technique  to  achieve  a  cylinder  40  microns  in 
diameter  and  20  microns  thick.  The  results  for  fused  silica  are  shown  in  figure  24. 

This  change  in  volume  can  be  converted  into  a  change  in  density  and  plotted  upon  previous 
measurements  for  fused  silica  (13).  Figure  25  provides  a  compilation  of  previous  studies  on  the 
density  changes  with  pressure  for  previous  methods.  The  diamonds  plots  represent  densified 
glass;  glass  that  has  been  pressurized  to  10  GPa  and  then  released  to  ambient  pressure. 
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Figure  24.  DAC  results  for  fused  silica  compared  to  previous  DAC  results 

published  in  (12)  (our  results  are  in  red).  Phase  lines  for  the  different 
quart/  structures  are  shown.  At  approximately  27  GPa,  the  sample  is 
thought  to  be  bridged  by  the  diamonds.  This  accounts  for  the 
deviation  seen  at  27-30  GPa. 


Figure  25.  Plot  of  density  changes  with  pressure  for  fused  silica.  Original  figure 
from  Wakabayashi  et  al.  (13). 
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Shown  in  figure  26  is  an  example  of  the  comparison  between  the  sample  at  4  and  33  GPa  for 
fused  silica.  The  gray  circle  on  the  4  GPa  side  corresponds  to  the  sample  diameter  at  33  GPa  and 
is  shown  for  comparison  purposes.  As  the  pressure  is  increased,  the  sample  hole  closes  due  to 
the  plastic  flow  of  the  gasket  material. 

In  addition  to  the  DAC  experiments,  neutron  diffraction  was  investigated  on  the 
nanoscale-ordered  materials  diffractometer  (NOMAD)  beamline  at  the  spallation  neutron  source 
(SNS)  at  Oak  Ridge  National  Laboratory.  This  technique  results  in  a  total  scattering  cross 
section,  which  is  integrated  over  the  entire  sample.  The  scattering  can  be  converted  into  a 
momentum  transfer  and  provide  a  “first  sharp  diffraction  peak”  (FSDP),  seen  in  figure  27.  The 
simplest  way  to  interpret  Q  is  to  realize  that  the  lower  the  Q  the  larger  the  interaction.  If  one 
thinks  about  the  network  structure  of  glasses,  the  scattering  can  be  broken  up  into  several  regions, 
namely,  the  SRO,  IRO,  and  long-range  order.  Glasses  have  rigid  SRO  in  the  form  of  a  silicon 
dioxide  tetrahedron  and  disorder  with  regard  to  the  random  nature  of  the  long-range  order.  IRO 
can  have  order  in  the  form  of  rings,  cages,  and  chains.  In  order  to  extract  the  interaction 
distances,  it  is  important  to  use  a  high  energy  source;  due  to  the  nature  of  the  Fourier  transform, 
the  higher  the  energy  the  better  the  resolution  on  the  low  Q  region. 


Figure  26.  Comparison  of  sample  upon  compression  in  the  DAC.  The  3.9-GPa 

sample  has  a  gray  sphere  to  indicate  the  area  of  the  sample  at  33.1  GPa 
to  assist  the  eye  in  comparison. 
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Figure  27.  Neutron  diffraction  spectra  of  fused  silica,  Borofloat,  and  SL  glasses. 
Intensities  are  not  corrected  for  the  boron  neutron  absorption  cross 
section.  The  location  of  the  first  sharp  diffraction  peak  is  indicated; 
shifts  to  higher  Q  correspond  to  decreases  in  intermolecular  interaction 
distances. 


To  further  study  the  effects  of  composition  on  pressure  and  ballistic  response,  samples  of  glass 
were  prepared  with  known  concentrations  of  secondary  molecules.  These  samples  were 
subjected  to  initial  study  with  the  laser  Raman  system  to  discern  if  there  were  noticeable 
differences  in  the  spectra  shown  in  figure  28.  The  spectra  were  obtained  using  a  532-nm  diode 
laser  (average  power  300  mW)  and  a  custom  Raman  spectrometer.  Preliminary  analysis  of  the 
results  seems  to  indicate  a  change  in  the  300-400  wavenumber  regions  with  decreasing  intensity 
as  the  composition  changes  with  each  sample.  The  spectra  were  normalized  to  a  higher 
wavenumber  region  (<  2500  cm^1).  The  samples  for  use  in  the  DAC  are  being  prepared  using 
the  FIB  so  similar  density  pressure  analysis  can  be  attempted. 
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Figure  28.  Raman  comparison  of  four  prepared  glass  samples  by  Dr.  Parimal 

Patel.  Peak  at  zero  is  the  laser  line,  partially  blocked  by  a  notch  filter 
that  extends  to  approximately  150  cm-1. 


4.2.10  Nanoindentation  Studies  of  Fused  Silica 

Nanoindentation  experiments  were  conducted  at  ARL  on  fused  silica  specimens  using  an  MTS 
Nanoindenter  XP  operated  in  continuous  stiffness  measurement  (CSM)  mode.  A  spherical 
indenter  with  a  radius  of  3  /im  was  used  to  indent  to  depths  approaching  2  //m  at  a  constant  strain 
rate  of  0.05/s.  The  hardness  values,  H,  were  calculated  from  the  maximum  loads,  Pmax,  and  the 
contact  area,  A,  at  the  maximum  indentation  depth  where  H  =  The  elastic  modulus  values 

were  calculated  using  the  Oliver  and  Pharr  method  (38)  from  the  measured  unloading  stiffness,  S 
(which  is  equivalent  to  the  slope  of  the  initial  unloading  curve),  as  follows: 


S  = 


(5) 


where  Eeff  is  a  function  of  the  Poisson’s  ratio  and  elastic  modulus  for  the  indenter  (v%,  EA  and 
material  of  interest  (v,  E)  defined  as 
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(6) 


1  _  1  -  v2  1  -  v‘l 

Eeff  E  +  Ei 

Following  the  indentation  experiments,  the  residual  indents  were  examined  with  a  NanoSEM  600 
SEM  operated  in  low-vacuum  mode  (which  is  used  to  accommodate  non-conductive  specimens.) 

Figure  29  shows  the  typical  load-displacement  curves  for  maximum  displacements  ranging  from 
500  to  2000  nm.  Figure  30a  plots  the  measured  hardness  as  a  function  of  indentation  depth  and 
indicates  there  is  some  indentation  size  effect  over  the  range  considered.  The  standard  deviation 
is  largest  for  the  smallest  indentation  depths  (500  nm).  The  elastic  moduli  decrease  with 
increasing  maximum  indentation  depth  (figure  30b);  however,  the  error  in  measurement  does  not 
follow  a  trend  with  depth.  The  SEM  examination  gives  insight  into  the  indentation  size  effect. 

We  are  unable  to  resolve  a  residual  impression  for  the  specimens  indented  to  a  depth  limit  of  500 
nm,  which  indicates  the  response  is  mostly  elastic  at  small  depths.  Figures  31-33  show  SEM 
micrographs  for  specimens  indented  to  depths  of  1000,  1500,  and  2000  nm,  respectively.  There 
is  a  noticeable  residual  indent  in  figure  3 1 ;  however,  there  is  evidence  of  fracture.  At  greater 
indentation  depths,  radial  cracks  (figures  32a  and  33)  and  classic  cone  cracks  (figure  32b)  are 
visible.  Further  investigation  is  required,  but  the  lower  hardness  and  modulus  values  measured  at 
greater  indentation  depths  could  result  from  indentation  cracks.  In  some  brittle  material  systems, 
“pop-ins”  or  discrete  jumps  in  displacement  during  indentation  are  found  in  the  load-displacement 
curves.  However,  no  pop-ins  are  observed  in  this  series  of  nanoindentation  tests. 
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Figure  29.  Load  vs.  displacement  in  fused  silica. 


13 


12.5 

12 


1 


11.5 

11 


* 


10.5 

10 

0  500  1000  1500  2000  250( 


Indentation  Depth  (nm) 


78 


76 


70 


68 


0 


500  1000  1500  2000 

Indentation  Depth  (nm) 


(a) 


(b) 


Figure  30.  (a)  Hardness  and  (b)  modulus  variations  in  fused  silica. 
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Figure  31.  SEM  micrographs  of  fused  silica  indented  to  1000  nm. 


Figure  32.  SEM  micrographs  of  fused  silica  indented  to  1500  nm  (a)  radial  cracks 
and  (b)  cone  cracks. 
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Figure  33.  SEM  micrographs  of  fused  silica  indented  to  2000  nm. 
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5.  Pairwise  Functional-free  Silica  Potential  From  First-principles 
Molecular  Dynamics  Simulation 


We  have  developed  a  new  short-range  pairwise  numerical  potential  for  silica,  as  described  in 
Izvekov  and  Rice  (14).  Here  we  present  the  key  features  of  the  model.  The  potential  is  derived 
from  a  single  AIMD  simulation  of  molten  silica  using  the  force-matching  (FM)  method,  with  the 
forces  being  represented  numerically  by  piecewise  functions  (splines)  (39).  The  AIMD 
simulation  is  performed  using  the  Bom-Oppenheimer  method  with  GGA  (BLYP)  and  XC 
functional.  The  new  effective  potential  shown  in  figure  34  includes  a  soft  repulsive  shoulder  to 
describe  the  interactions  of  oxygen  ions.  The  new  potential,  despite  being  short-ranged  and 
derived  from  single-phase  data,  exhibits  a  good  transferability  to  silica  crystalline  polymorphs 
and  amorphous  silica.  The  importance  of  the  0  —  0  soft  repulsive  shoulder  interaction  on  glass 
densification  under  cold  and  shock  compressions  is  assessed  from  MD  simulations  of  silica  glass 
under  room  and  shock  Hugoniot  conditions,  respectively,  and  shown  in  figure  35.  Results  from 
these  simulations  indicate  that  the  appearance  of  oxygen  complexes  (primarily  pairs)  occurs  at 
8-10  GPa,  and  under  cold  compression  conditions  becomes  notable  at  40  GPa,  essentially 
coinciding  with  the  transition  to  a  Si  sixfold  coordination  state.  An  analysis  of  changes  in  system 
structure  in  compressed  and  shocked  states  reveals  that  the  O  ions  interacting  through  the  soft 
repulsive  shoulder  potential  in  denser  states  of  silica  glass  may  create  a  mechanical  multi- stability 
under  elevated  pressures,  and  thus  contribute  to  the  observed  anomalous  densification. 

At  pressures  P  <  8-10  GPa,  the  densification  by  the  FM  model  occurs  predominantly  due  to 
structure  deformation  and  topological  reorganization  rather  than  a  weakening  of  the  0  —  0 
repulsion  due  to  the  soft  repulsive  shoulder  of  the  FM  potential.  The  observed  disagreement  of 
the  FM  cold  and  Hugoniot  density-pressure  curves  at  low  (P  =  0-10  GPa)  pressures  with  the 
experiment  data  (figure  35)  can  be  explained  by  the  deficiency  of  the  FM  model,  which  is 
pairwise  and  central,  regarding  angle-bending  terms.  The  extension  of  the  pairwise  FM  model 
with  angle-bending  forces,  which  are  three-body  forces,  is  necessary  to  improve  performance  of 
the  FM  model  at  low  pressures  in  which  densification  is  driven  by  mostly  structural  and 
topological  changes.  This  work  is  in  progress. 
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Figure  34.  Effective  atom-atom  forces  [panel  (a)]  and  corresponding  potentials 
[panel  (b)]  in  liquid  Si()2  generated  through  the  FM  method  as 
functions  of  interatomic  separation:  0  —  0  (black),  Si  —  O  (red),  and 
Si  —  Si  (green).  The  dashed  line  corresponds  to  the  model  without 
the  soft  repulsive  shoulder  (FM-ns  model).  In  panel  (b),  the  dotted 
line  corresponds  to  the  FM  model  with  a  weaker  repulsive  shoulder 
(FM-ws  model)  and  the  dot-dot-dashed  lines  indicate  the  variation  of 
the  O  —  O  repulsion  in  the  FM  procedure  with  a  length  of  ab  initio 
reference  trajectories  as  discussed  in  {14). 
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Figure  35.  Density  vs.  pressure  at  T  =  298  K  (solid  lines)  and  along  the 
Hugoniot  locus  (dashed  lines)  from  simulations  of  glass  sample 
structure  using  the  FM  (cyan),  FM-ns  (green)  [panel  (a)],  Pedone 
(blue)  [panel  (b)],  and  FM-ws  (magenta)  [panel  (c)]  models. 
Experimental  EOS  obtained  by  cold  compression  (black  squares)  and 
by  shock  compression  (red  circles)  are  from  (15)  and  (16), 
respectively.  Insert  to  panel  (c)  compares  the  298  K  (solid)  and 
Hugoniot  (dashed)  EOS  from  simulations  of  the  ambient  density 
(FM-1)  and  densified  (FM-h)  glass  samples  using  the  FM  models. 
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5.1  Molecular  Dynamics  Modeling  of  Glass  Nanoindentation 


MD  methods  have  been  used  to  study  nanoindentation  for  a  number  of  material  systems,  i.e., 
metals  (40),  ceramics  (41),  glasses  (42),  and  energetic  materials  (43).  Length  scales  for  MD 
simulations  can  be  made  comparable  to  those  in  the  experiments  of  Nomura  et  al.  (42),  although 
the  time  scales  and  strain  rates  differ  significantly.  The  advantage  of  MD  simulations  over 
nanoindentation  experiments  is  the  capability  to  provide  atomistic  detail  of  numerous  properties, 
including  stress  distribution  and  structural  information. 

We  have  previously  reported  on  MD  simulations  of  nanoindentation  performed  on  a  large-scale 
fused  silica  system  (44).  Here  we  extend  the  analysis,  including  a  calculation  of  the  hardness, 
and  compare  to  experimental  results.  As  described  in  Gazonas  et  al.  (44),  all  the  simulations 
were  performed  using  the  Large-scale  Atomic/Molecular  Massively  Parallel  Simulator 
(LAMMPS)  (45),  with  the  pairwise  potential  recently  developed  using  FM  techniques  described 
in  section  5.  The  simulation  cell  is  29.9  x  29.9  x  17.8  nm,  containing  1,160,952  atoms  (386,984 
a  —  Si02  molecules).  Once  the  system  is  equilibrated,  following  the  annealing  schedule 
described  in  Pedone  et  al.  (46),  a  spherical  indenter  with  a  radius  of  9  nm  is  introduced  in  the 
z-direction.  The  indenter  interacts  with  atoms  in  the  simulation  cell  via  a  force  of  magnitude, 


F  =  —K(r  -  Rf 


(7) 


where  K  is  a  force  constant,  R  is  the  radius  of  the  indenter,  and  r  is  the  distance  from  the  atom  to 
the  center  of  the  indenter.  Periodic  boundary  conditions  are  implemented  in  the  x-  and  y- 
directions;  the  indented  surface  remains  free,  and  approximately  33,000  atoms  on  the  opposite 
surface  are  held  fixed  to  ensure  that  the  system  remains  stationary  during  indentation.  The 
simulation  is  performed  in  the  microcanonical  ensemble  (NVE),  with  a  timestep  of  2.0  fs.  It  has 
recently  been  determined  that  this  timestep  is  too  large  for  simulations  using  the  FM  potential. 
The  dynamics  of  the  system  are  not  properly  captured  when  using  a  large  integration  timestep, 
and  the  differences  between  our  results  and  those  from  experiment  are  due  partly  to  this  large 
timestep,  as  is  discussed. 

Figure  36  shows  the  complete  loading  and  unloading  curve,  as  well  as  the  cross  sections  of  the 
corresponding  atomic  configurations  at  the  maximum  load  (figure  36b),  and  after  complete 
unloading  (figure  36c).  The  hardness,  H,  is  defined  as, 
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Figure  36.  (a)  Force  vs.  indenter  depth  for  a  complete  loading  and  unloading 

cycle.  Atomic  configuration  (b)  at  the  point  of  maximum  loading,  and 
(c)  after  complete  unloading. 


P 

TT  ma%  / Q\ 

H  =  —  •  (8) 

where  Pmax  is  the  maximum  applied  load  and  A  is  the  contact  area.  The  maximum  load  is 
determined  from  the  loading  curve  in  figure  36a  by  fitting  the  latter  portion  of  the  data  to  a  line 
and  determining  the  load  corresponding  to  the  maximum  depth.  The  value  we  calculate  is  1.10  n 
N.  To  calculate  the  contact  area,  we  use  a  method  similar  to  that  described  in  Chen  and  Ke  (47). 
First,  the  atoms  in  contact  with  the  indenter  are  identified  by  calculating  the  distance  between 
each  atom  and  the  center  of  the  indenting  sphere.  Those  atoms  appear  in  red  in  figure  37,  along 
with  their  projection  onto  the  xy-plane.  The  area  of  a  circle  in  the  xy-plane  that  contains  all  the 
projections  is  defined  as  the  contact  area.  We  calculate  an  area  of  14900  A3,  resulting  in  a 
hardness  value  of  7.38  GPa.  This  differs  from  the  value  we  previously  reported  (44)  due  to  the 
more  refined  method  of  calculating  the  contact  area.  Furthermore,  it  is  not  in  agreement  with  the 
experimental  value  of  10  GPa  reported  in  Miyake  et  al.  (48),  nor  with  the  depth-dependent 
hardness  values  determined  from  our  own  nanoindentation  experiments  (see  figure  30a).  As 
mentioned  above,  the  disagreement  could  be  attributable  to  the  too  large  timestep,  but  the  effect  of 
sample  size  should  also  be  explored  to  ensure  the  size  of  the  sample  is  not  adversely  affecting  the 
hardness  calculation. 
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Figure  37.  The  atoms  that  are  in  contact  with  the  indenter  are  shown  in  red,  and 
their  projection  onto  the  xy-plane  is  shown  in  black. 


Nanoindentation  experiments  have  shown  that  the  unloading  curve  can  be  approximated  by  the 
power  law  relation  (38), 


P  =  a(h  —  hf)m 


(9) 


where  a  and  m  are  fitting  constants,  and  the  final  depth,  hf,  is  the  permanent  depth  of  penetration 
after  the  indenter  is  fully  unloaded  (see  figure  36).  We  have  fit  the  unloading  curve  in  figure  36a 
to  this  equation,  with  hf  included  as  a  fitting  parameter,  resulting  in  a  =  0.9720  //N/nnT",  m  = 
1.445,  and  hf  =  2.889  nm.  The  value  of  hf  is  in  agreement  with  our  data,  and  the  value  of  m  for 
a  spherical  indenter  is  expected  to  be  1.5  (49),  close  to  the  value  we  obtain. 

The  value  of  a  differs  significantly  from  the  experimental  value  of  50.0  //N/nm'"  for  fused 
silica  (38),  although  that  was  obtained  with  a  Berkovich  indenter.  The  discrepancy  in  a  is  not 
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surprising  as  the  loading  depths  reached  in  nanoindentation  experiments  are  much  larger  than  in 
our  simulation,  and  the  strain  rates  are  significantly  lower  than  can  be  achieved  in  MD 
simulations.  However  the  shape  of  the  curve,  indicated  by  the  value  of  m,  is  consistent  with 
experimental  values,  despite  the  timestep  being  too  large.  In  addition  to  repeating  the  above 
simulation  with  a  smaller  timestep,  we  plan  to  perform  nanoindentation  on  larger  systems  to  study 
size  effects  on  the  mechanical  response  of  fused  silica.  We  will  also  perform  nanoindentation 
simulations  using  a  flat  punch  indenter  (for  which  m-  1  in  equation  9)  to  directly  compare  with 
the  dimensional  analysis  results  presented  in  section  6.  Finally,  the  force  match  potential  will  be 
extended  to  include  three-body  forces  as  well  as  additional  atom  types  in  order  to  model 
borosilicate  glasses. 


6.  Dimensional  Analysis  of  Hertzian  Cone  Crack  Development  in  Brittle 
Elastic  Solids 


In  this  section,  we  outline  the  methodology  for  the  derivation  of  similarity  relationships  for  elastic 
solids  exhibiting  stable  crack  growth  induced  by  rigid  indenters;  as  will  be  shown  in  a  subsequent 
section,  such  relationships  are  invaluable  for  verification  of  computational  models  that  involve 
crack  propagation.  The  application  of  dimensional  analysis  to  the  development  of  stable  cone 
cracks  appears  to  have  been  first  addressed  by  Roesler  (77)  and  Benbow  (18),  and  later  more 
fully  analyzed  by  Barenblatt  (50).  The  universal  dimensionless  relationship: 


■  <io) 

relates  the  width  D  of  the  base  of  a  cone  crack  induced  by  a  flat  punch  (cylindrical  indenter)  to 
the  load  P,  and  cohesive  modulus  K  of  the  medium  (figure  38).  Cone  crack  development  in 
glass  (figure  39)  confirms  the  universal  scaling  law  (equation  10)  illustrated  in  the  log-log  plot  of 
figure  40.  In  addition  to  the  relevant  physical  parameter  illustrated  in  figure  38,  Poisson’s  ratio, 
v,  and  the  modulus  of  cohesion  K  play  an  important  role  in  describing  the  mechanics  of  the 
indentation  problem.  According  to  similarity  methods  described  by  Sedov  (51)  and  others,  since 
there  are  five  relevant  physical  variables,  and  three  fundamental  dimensions  of  M,  L,  and  T,  this 
results  in  two  possible  dimensionless  products,  which  can  be  formed  from  the  physical  variables. 
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Figure  38.  Dimensional  analysis  for  the  axisymmetric  Hertzian  cone  crack 
problem  . 


Figure  39.  Hertzian  cone  crack  in  SL  glass  induced  by  a  flat  punch  (cylindrical 
indenter)  after  (17)  . 
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Figure  40.  Hertzian  cone  crack  data  for  fused  silica,  and  experimental 
confirmation  of  scaling  law  from  Benbow  (18). 


The  two  dimensionless  groups  can  be  derived  for  this  problem  by  first  writing  the  relevant 
physical  parameters  in  the  following  multiplicative  form: 


(. D)kl(d0)k2(P)k3(K)k*(v)ks  =  1 


and  the  expressing  equation  1 1  in  terms  of  the  fundamental  dimensions: 

(L)kl(L)k2(MLT-2)k3(ML-1/2T~2)k4{  l)*5  =  1 


The  fundamental  dimension  terms  in  equation  12  can  be  factored  as  follows: 


(M)k3+ki(L)kl+k2+k3-ki/2(T)-2k3-2k4  =  i 


resulting  in  the  following  system  of  equations: 


(ID 


(12) 


(13) 


47 


5~  —  0 

Aq  “I-  /^2  “1“  /?3  —  /u4/2  =  0 

— 2k3  —  2A;4  =  0 


(14) 


that  are  solvable  using  the  constants  that  appear  in  table  11. 


Table  11.  Constants  appealing  in 
dimensional  analysis  . 


h 

&2 

h 

h 

7Ti 

1 

0 

-2/3 

2/3 

1 

VT2 

0 

1 

-2/3 

2/3 

1 

The  dimensionless  groups,  also  known  as  7r-groups  from  the  Buckingham  7r  theorem  (52)  can  be 
written  as 


7T 1  =  DP~2/3K'2/3U  , 

7T2  =  doP~2/3K2/3v  .  (15) 

Since  the  7r-groups  in  equation  15  are  dimensionless,  they  can  be  equated  and  solved  for  the  width 
D  of  the  base  of  the  cone  crack, 

/  p  \  2/3  js 

d=\k)  "/w>(;p)2/3’!/)  '  <16) 

Thus,  a  log-log  plot  of  D  versus  P,  for  D  »  d0,  will  have  a  slope  of  2/3,  as  it  appears  that  the 
contribution  from  the  third  function  in  equation  16  is  negligible  (see  figure  40).  The  fracture 
scaling  law  can  also  be  used  to  validate  finite  element,  peridynamic,  material  point,  MD,  and 
other  computational  methods  currently  in  use  for  simulating  fracture  resulting  from  indentation 
experiments. 

For  fractures  that  develop  under  two-dimensional  (2-D)  states  of  stress  such  as  the  problem  of  the 
symmetric  wedging  of  a  thick  plate  (79),  the  wedging  force  can  be  approximated  by  two  equal 
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and  opposite  concentrated  tractions  P  per  unit  length  of  plate  thickness  as  illustrated  in  figure  41. 
Using  the  dimensional  analysis  methods  outlined  earlier  in  this  section,  a  dimensionless  relation 
can  be  derived,  for  2-D  stress  and  deformation  fields,  which  only  depends  on  crack  length  l,  load 
P,  and  cohesive  modulus  K  (19): 


Figure  4 1 .  Idealized  2-D  problem  of  wedge-induced  crack  propagation  in  a  thick 
plate  after  Barenblatt  (19). 


This  result  is  used  to  validate  our  peridynamic  code  simulations  of  fracture  in  section  10.1. 


7.  The  Mechanics  of  Indentation 


In  this  section,  we  outline  the  exact  solution  to  the  axisymmetric  indentation  problem  into  an 
elastic  halfspace,  otherwise  known  as  the  Boussinesq  problem,  that  is  relevant  to  the  dimensional 
analysis  described  in  section  6.  The  exact  elastic  solutions  provide  a  means  to  quantitatively 
validate  computational  simulations,  prior  to  the  development  of  permanent  densification  and 
fracture  of  the  glasses  under  consideration  in  this  research.  Boussinesq  (53)  and  Love  (54)  were 
among  the  first  to  consider  such  a  problem,  but  Sneddon  solved  the  problem  for  the  flat-ended 
punch  (55,  56)  in  cylindrical  polar  coordinates  using  Hankel  transforms  that  resulted  in  the 
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solution  of  dual  integral  equations  (57).  The  solutions  to  the  field  equations  are  listed  below  as 
they  appear  in  Sneddon’s  (25)  textbook, 


CTg  = 


(7  r 


Ur  = 


Ur  = 


^  (4  -  ^Wi) 

7t(A  +  2  fi) 

2‘  +  JS) 


7 r 


Or  = 


4/ie  (C^2  +  Jj)  (A  +  p) 
7ra(A  +  2/i) 
A<^Jlp,e(\  T  /i) 


7ra(A  +  2/i) 

4CJf  iljh  (J»  - 


7ra(A  +  2/r)  nap(\  +  2/i) 

4/ie  ( (2A  +  /i)  —  C^2  +  A4)) 


00  = 


7ra(A  +  2/i) 


(18) 


In  these  equations,  the  dimensionless  radial  coordinate  p  —  r/a  is  the  physical  radial  distance  r 
normalized  by  the  indenter  radius  a,  and  the  dimensionless  depth  coordinate  (  =  z/ a  is  the 
physical  depth  z  normalized  by  the  indenter  radius  a.  e  is  the  indentation  depth,  and  the  Lame 
parameters  are  given  by  A  and  /j.  The  radial  stress  ar  is  not  given  in  Sneddon’s  text  but  can  be 
derived  readily  by  subtraction  of  og  from  ar  +  og  from  equations  18  resulting  in 


(J  r 


4 :/te(  Jg/t 


(A  +  /X)(p(j0-CJ20)  +  C./i)) 
nap(\  +  2/i) 


where  Bessel  functions  J™  are  defined  as 


(19) 


jm  _ 

Jn  — 


pn  1e  sin(p)  Jm(jpp)  dp 


or  as  the  imaginary  part  of  the  integral 


(20) 
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(21) 


J?  = 


pn-le-p(C -i)jm(pp)dp 


We  derive  an  explicit  representation  for  equation  21, 


K  =  a  +  n)  2fi  i(m  +  n  +  1); m  +  1;  )  ,(22) 

where  (p,  (')  >  0, 5ft(m  +  n)  >  0,  S  is  the  imaginary  part  of  the  term  in  brackets,  m,  n  G  Z,  and 
2Fi  is  the  regularized  hypergeometric  function  2.F\(a,  b]  c;  z) /r(c). 

All  solutions  listed  in  equation  18  are  exact  and  expressible  in  closed-form  except  for  uz,  which 
requires  the  evaluation  of  J®.  The  solution  for  uz,  which  involves  Jq  is  written  in  terms  of 
T(m  +  n)  =  (m  +  n  —  1)!,  which  is  an  undefined  quantity  for  m  —  0,  n  —  0  as  this  violates  the 
condition  5ft  (m  +  n)  >  0  in  equation  22  and  probably  explains  why  the  full  plane  solution  for  uz 
is  not  provided  in  any  publications  we  are  aware  of  (25),  (55),  (5(5),  (58)*.  Sneddon’s  solutions 
can  be  rewritten  more  succinctly  using  the  substitution,  A  =  and  by  normalizing  the  stress 
components  by  the  mean  pressure  pm  =  ^a+v]  ,  and  the  displacement  components  by  the 
indentation  depth  e,  resulting  in 


(7  r 
Pm 


ur  _  (~C^1  +  (1  ~  2^)Jq) 

e  7r(l  —  v) 

uz  (2(l-z/)J0°  +  CJ1°) 

e  7r(l  —  v) 

Trz  (J2 

Pm  2 

—L  —  _1  />  7O  |  t0\ 

-  2(W2+-A) 

ym  ^ 

^  =  -  Ac  (topj*  -  jd + (1  -  2  v)4\ 
=  ~l -  C (t>4  -  J l)  -  (1  -  20^1 


(23) 


Note  that  the  analytical  solutions  in  equation  23  for  elastic  displacements  and  stresses  in  the 

*The  first  author  became  aware  of  an  explicit  expression  for  Jq  in  appendix  2  of  Barquins  and  Maugis  (59)  taken 
from  Dahan  (60),  prior  to  the  completion  of  this  final  report. 
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halfspace  are  functions  of  Poisson’s  ratio  v,  which  substantiates  the  choice  of  v  alone  to  describe 
the  elastic  properties  of  the  halfspace  in  the  dimensional  analysis  described  in  section  6. 

7.1  Principal  Stresses  in  a  Halfspace  Under  Axisymmetric  Indentation 

The  principal  stresses  G\  >  a2  >  a3  in  a  halfspace  subjected  to  indentation  by  an  axisymmetric 
cylindrical  indenter  (figures  42-  44)  can  be  determined  by  finding  the  eigenvalues  of  the 
following  matrix: 


(  ar  0 

0  (To 
\  ~i"rz  0 


T~rz  ^ 
0 

) 


which  are 


(24) 


°1  =  \  +  CTz  +  V7 ar  -  2oraz  +  4rr2,  +  ,  (25) 

°"3  =  2  (a'r  +  CTz~  “  2(jr(jz  +  4rr2*  +  afj  . 

We  also  plot  the  normalized  vertical  displacement  —  of  the  halfspace  (figure  45)  using  both  the 
exact  solution  from  Sneddon’s  text  (25)  as  well  as  an  approximate  solution  derived  by  using  both 
a  series  expansion  approximation  of  the  Bessel  function  J0  (pp)  near  the  origin  and  the  asymptotic 
expansion  of  J0  (pp)  for  large  p  that  appears  in  equation  1 8  or  equation  23  in  the  term  involving 


Jo=  e  Cpsinc (p)J0(pp)dp  .  (26) 

Jo 

We  note  that  Sneddon’s  text  (25)  and  other  references  on  this  topic  provide  the  solution  only  for 
the  halfspace  surface  displacement  field,  because  of  the  difficulty  in  evaluation  of  Jq  ;  the  topic  of 
determining  the  displacement  field  —  throughout  the  entire  halfspace  will  be  the  subject  of  an 
upcoming  publication.  Finally,  we  use  the  exact  indentation  solution  for  the  surface  displacement 
of  a  halfspace  found  in  Sneddon’s  text  (25)  (also  plotted  in  figure  45)  to  validate  our 
axisymmetric  peridynamic  code  in  section  10.3. 
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Figure  42.  ^p-  contours  in  the  halfspace  for  =  0-1679,  where 

A  =  15.833  GPa,  //  =  31.3  GPa  using  the  Lame  parameters  for  fused 
silica  derived  from  Scheidler  (20). 
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Figure  43.  ^p-  contours  in  the  halfspace  for  v  =  =  0.1679,  where 

A  =  15.833  GPa,  fj,  =  31.3  GPa  using  the  Lame  parameters  for  fused 
silica  derived  from  Scheidler  (20)  . 
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Figure  44.  j3-  contours  in  the  halfspace  for  =  2(\+IJ,)  =  0-1679,  where 

A  =  15.833  GPa,  fj,  =  31.3  GPa  using  the  Lame  parameters  for  fused 
silica  derived  from  Scheidler  (20)  . 
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Figure  45.  Deformed  surface  of  the  halfspace  for  u  =  2(\+^)  =  0.1679,  where 
A  =  15.833,  /r  =  31.3  using  the  Lame  parameters  for  fused  silica 
derived  from  Scheidler  (20)  . 
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8.  Quantum  Mechanics  Modeling  of  Densification  and  Bulk  Modulus  of 
Silica  Versus  Pressure 


A  detailed  understanding  of  the  densification  process  and  structural  changes  in  amorphous  solids 
under  pressure  is  appealing  for  both  experimental  and  simulation  work.  One  of  the  interesting 
questions  is  about  the  nature  of  the  structural  transformation  between  low  and  high  density 
amorphous  phases.  To  model  the  structure  under  pressure  from  first  principles,  we  used  models 
with  different  densities,  number  of  atoms,  and  different  ring  statistics.  We  used  two  different 
methods  to  generate  the  random  connected  networks  (1)  a  Monte-Carlo  bond  switching  model 
(72-  and  1 14-atom  models)  and  (2)  MD  simulated  annealing  of  melted  silica.  The  ring  statistics 
could  be  described  by  plots  of  ring  size  distribution,  as  seen  in  figure  46.  The  faster  the 
quenching  of  the  melt  is,  the  wider  is  the  distribution  of  the  ring  sizes.  Here  is  an  example  of 
slow  quenching  with  ring  distributions  from  4  to  8  member  rings,  which  corresponds  to  the 
quartz-like  structure  of  the  1 14-atom  model.  Two  other  models  with  72  and  192  atoms  have 
wider  distribution  of  the  ring  sizes  from  3  to  12  member  rings;  for  reference,  quartz  has  only  6 
member  rings.  Three  to  four  member  rings  have  a  low  concentration,  but  play  an  important  role 
because  they  correspond  to  the  most  reactive  sites.  The  angles  between  the  Si  —  O  —  Si  atom, 
and  the  O  —  Si  —  O  atom  distributions  also  convey  important  information  about  structural 
changes  under  certain  pressure  (figure  47). 
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Figure  46.  Ring  distribution  for  the  1 14-atom  model. 
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Figure  47.  Two  types  of  angle  distribution  in  the  1 14-atom  model. 


By  relaxing  the  internal  coordinates  under  compressive  or  tensile  pressures  at  normal  conditions, 
we  found  that  the  Si  —  O  —  Si  angle  ranges  between  130°  and  180°,  while  the  tetrahedral  units 
are  preserved.  This  result  is  in  good  agreement  with  previous  theoretical  and  experimental 
results  of  Mauri  et  al.  (61).  Clear  signs  of  structural  transformations  in  silica  under  the  pressure 
may  be  seen  from  the  calculated  mass  density  (figure  48).  We  did  the  calculations  of  density 
after  optimizing  the  shape  and  volume  of  the  unit  cell  without  projections  in  real  space.  The 
analysis  used  the  Projector  Augmented  Wave  (PAW)  method  implemented  in  the  Vienna  ab  initio 
simulation  package  (VASP)  code  described  in  Kresse  and  Hafner  (62). 

In  all  three  models,  one  may  see  two  slopes  of  density,  which  might  be  related  with  different  types 
of  structural  changes.  One  occurs  up  to  20-30  GPa,  and  another  region  is  from  30-50  GPa. 

Some  signatures  of  the  two  phases  may  be  seen  from  the  x-ray  absorption  experiments  of  Sato 
and  Funamori  (15)  (figure  49). 

To  understand  what  is  so  special  with  these  two  stages  of  structural  transformation  in  silica,  we 
did  angle  distribution  analysis  under  pressures  similar  to  that  depicted  in  figure  47.  It  turned  out 
that  in  the  first  region  there  is  not  much  change  in  O  —  Si  —  O  distribution,  but  discernible 
changes  in  Si  —  O  —  Si  distribution  indicating  that  up  to  20-30  GPa  there  is  a  squeezing  of  space 
between  tetrahedra  and  not  much  distortion  of  tetrahedrons.  Each  O  atom  has  two  nearest 
neighbors;  each  Si  atom  has  four  nearest  neighbors.  At  pressures  higher  than  30  GPa,  silica 
becomes  very  dense,  and  both  the  Si  —  O  —  Si  and  O  —  Si  —  O  angles  change,  revealing 
transformation  in  both  the  tetrahedra  and  the  space  between  them. 
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Figure  48.  Density  of  fused  silica  as  a  function  of  pressure  for  72-atom  (T), 
114-atom  (▼),  and  192-atom  (Blmodels. 


Figure  49.  Experimental  density  of  fused  silica  after  Sato  and  Funamori  (75). 
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At  40-50  GPa,  severe  distortion  of  the  tetrahedrons  resulting  in  formation  of  sixfold-coordinated 
Si  atoms  (figure  50).  The  density  functional  theory  (DFT)  observation  of  the  sixfold-coordinated 
atoms  confirms  the  assumption  suggested  in  Sato  and  Funamori  (75). 

A  manifestation  of  the  two  phases  of  silica  densification  might  exhibit  the  unusual  behavior  of  the 
elastic  constants,  particularly  the  bulk  modulus.  We  calculated  the  bulk  modulus  from  the 
stress-strain  relationship  for  1 14-atom  model.  The  bulk  modulus  of  fused  silica  generally 
increases  with  pressure,  but  unusual  behavior  of  the  pressure  dependence  up  to  20  GPa  may  be 
related  to  the  densification  of  the  space  between  the  tetrahedra  (figure  51).  A  second  region  after 
20  GPa  with  significant  increase  of  bulk  modulus  corresponds  to  the  densely  packed  tetrahedra. 


Figure  50.  Fused  silica  structure  under  50  GPa  with  sixfold-coordinated  Si  atoms. 


Figure  51.  Bulk  modulus  of  fused  silica  as  a  function  of  pressure. 
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8.1  Force  Matching  Pair  Potentials  for  Borosilicate  Glasses 


To  generate  pair  potentials  for  borosilicate  glasses,  we  used  the  FM  method  as  applied  by  Izvekov 
and  Rice  (14).  The  method  is  based  on  DFT  calculations  of  trajectories  in  the 
Born-Oppenheimer  approximation  at  5000  K.  The  pair  potentials  for  S  —  B  and  B  —  O  generated 
using  the  method  are  shown  in  the  figure  52. 

Pair  Si  —  Si,  Si  —  O,  and  0  —  0  potentials  are  similar  to  those  used  for  pure  silica  MD 
calculations.  The  numerical  pair  potentials  include  to  a  certain  extent  many  body  interactions 
since  they  were  generated  based  on  force  matching  of  DFT  calculations. 


(a)  (b) 

Figure  52.  Pair  potentials  for  (a)  Si  —  B  and  (b)  O  —  B. 


8.2  Simulation  of  Vibration  Spectra  of  Amorphous  Si02 

Given  the  large  variety  of  structural  building  blocks  of  amorphous  Si02,  there  is  s  need  for  a 
fundamental  description  of  the  random  lattice  vibrations,  particularly  those  that  are  IR  and  Raman 
active.  In  this  section,  IR  and  Raman  spectra  were  computed  using  DFT  for  both  crystalline  and 
amorphous  Si02.  In  an  insulating  crystal,  the  frequency  and  the  intensity  of  the  Raman  peaks 
are  determined  by  the  zone-center  phonon  frequencies  and  the  Raman  tensor.  The  phonon 
frequencies  are  determined  by  the  dynamical  matrix,  dielectric  constant,  and  Born  effective 
charges.  The  Bom  effective  charge  tensor  of  an  ion  is  the  partial  derivative  of  the  macroscopic 
polarization  with  respect  to  a  periodic  displacement  of  all  the  periodic  images  of  that  ion  at  zero 
macroscopic  electric  field.  Calculations  are  done  with  the  CASTEP  code  (63)  using 
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norm-conserving  pseudopotentials  and  plane  waves  with  600  eV  cutoff.  The  calculated  IR 
spectrum  reproduces  the  main  peaks  observed  in  a-quartz  as  shown  in  figures  53  and  54. 


Wavenumber  (cm1) 


Figure  53.  Calculated  IR  spectrum. 
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Figure  54.  IR  spectrum  was  registered  on  Nicolet  6700  FT-IR  spectrometer. 


The  atomistic  and  charge  relaxed  model  and  scaled  arrows  indicating  the  main  vibration  modes 
are  shown  figure  55. 
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Figure  55.  (a)  Relaxed  atomic  structure  and  charge  total  distribution  isosurface  at 
value  0.36  together  with  arrows  indicating  vibration  modes,  and  (b) 
distribution  of  local  self  consistent  potential. 


The  Raman  spectrum  calculated  in  linear  response  approximation  is  shown  in  figure  56.  The 
calculated  spectrum  reproduces  two  experimental  frequencies  at  466  and  205  cm-1.  The 
calculation  of  vibration  spectra  of  amorphous  Si02  depends  on  the  atomic  and  topological  model 
of  silica,  ring  size,  and  density.  The  72-atom  model  obtained  using  the  bond  switching 
algorithm  (64)  with  random  angles  and  no  broken  bonds  was  used  for  computing  the  vibration 
spectra.  The  structure  has  a  distribution  of  4-8  rings  (figure  57)  and  mass  density  of  2.145  g/cm3. 

The  calculated  IR  spectrum  of  the  silica  model  is  shown  in  figure  58.  The  main  vibration  modes 
corresponding  to  the  peaks  of  the  intensity  are  shown  by  scaled  arrows  in  figure  59a. 

It  is  evident  that  there  is  more  charge  within  random  rings  than  in  quartz  and  the  main  vibration 
modes  are  associated  with  atoms  located  in  larger  rings.  Peaks  around  474  and  1188  cm-1  are 
seen  in  simulated  curve  of  figure  53.  The  IR  band  at  800  cm-1  can  be  assigned  to  Si  —  O  —  Si 
symmetric  stretching  vibrations,  whereas  the  IR  band  at  474  cm-1  is  due  to  O  —  Si  —  O  bending 
vibrations.  The  calculated  Raman  spectrum  of  silica  is  shown  in  figure  60. 

In  the  experimental  Raman  spectra  of  silica  shown  in  figure  61,  three  regions  are  typically 
specified  (65).  The  broadband  around  490  cm^1  (Dl  band),  peak  around  600  cm-1  (D2  band) 
(attributed  to  symmetric  breathing  mode  of  3-4  member  rings),  and  peak  around  800  cm-1 
(assigned  to  the  Si  —  O  stretching). 

The  calculated  Raman  spectra  from  figure  60  reproduce  the  three  main  regions  of  the  spectrum  in 
figure  61  and  the  methodology  of  calculation  will  be  used  for  identifying  signatures  of 
non-standard  coordination  of  silica  atoms  at  high  pressures. 
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Figure  56.  Raman  intensity  of  quartz  (a)  calculated  and  (b)  experimental  (21). 
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Figure  57.  Distribution  of  rings  in  72-atom  model.  The  six  member  rings  have 
the  highest  concentration  in  the  model. 


64 


6000 


Wavenumber  (1/cm) 


Figure  58.  IR  spectrum  of  the  72-atom  silica  model. 


(a) 


(b) 


Figure  59.  (a)  The  main  vibration  modes  in  silica  are  shown  by  scaled  arrows. 
The  semi-transparent  clouds  represent  isosurface  of  total  charge 
distribution  at  a  value  of  0.35  and  (b)  experimental  Fourier  transform 
infrared  (FTIR)  spectrum  of  precipitated  amorphous  Si02  (22). 
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Figure  60.  Raman  spectra  of  72-atom  model  of  silica  using  10  cm-1  smearing,  10 
K  temperature,  and  Ar  laser  wavelength. 


Figure  61.  Experimental  Raman  spectra  of  silica  at  ambient  conditions  (lowest 
curve). 
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9.  Continuum  Equation-of-state  Model  Development 


Established  processes  to  develop  constitutive  models  for  use  in  large-scale  analysis  and  design 
simulations  are  heavily  influenced  by  the  vast  body  of  work  on  metals.  Typically,  a  model  is  split 
into  an  EOS,  relating  pressure,  density,  energy,  and  temperature,  and  a  strength  model,  which 
relates  the  deviatoric  part  of  the  stress  tensor  to  the  time  history  of  the  strain  rate,  temperature, 
and  state  variables.  This  is  the  context  of  the  glass  EOS  in  Gazonas  et  al.  (66),  where  a 
polyamorphic  model  was  used,  in  conjunction  with  high  and  low  pressure  EOS  models,  to 
represent  the  permanent  volume  change  in  glass  at  high  pressures  (in  excess  of  10  GPa). 

The  model  assumptions  are  evaluated  by  comparing  simulation  results  with  velocity  profiles 
obtained  from  gas  gun  experiments  by  Alexander  et  al.  (23).  The  experimental  results  are  shown 
for  three  different  impact  velocities  in  figure  62a.  The  initial  velocity  ramp  from  0  to  0.2  km/s 
results  from  the  decreasing  modulus  with  pressure,  which  is  not  captured  by  the  current  model. 

At  velocities  of  approximately  0.6  km/s,  there  is  a  kink  associated  with  the  polyamorphic 
transformation.  The  high  velocity  curve  (208  m/s  impact  velocity)  rises  abruptly  after  the 
transformation  point,  and  caps  at  a  velocity  consistent  with  the  initial  velocity  and  shock 
impedance  of  the  flyer.  The  transformation  is  nearly  over  driven  at  this  high  velocity,  in  that  there 
is  very  little  shift  between  the  two  vertical  sections  of  the  curve.  The  lower  curve  (137  m/s) 
shows  a  gradual  rise  associated  with  the  kinetics  of  the  transformation.  It  reaches  the  peak 
velocity  by  the  end  of  the  plot.  The  profile  from  the  intermediate  velocity  experiment  retains 
some  curvature  due  to  the  transformation  kinetics,  but  the  steady  velocity  is  attained  quickly. 

Results  from  using  the  initial  model,  with  separate  elastic-plastic  and  polyamorphic 
transformations,  in  gas  gun  simulations  yields  velocity  profiles  displaying  a  three-wave  structure 
(figure  62b)  rather  than  the  two  waves  seen  in  the  experiments.  The  model  kinetics  were 
assumed  to  be  fast  in  these  simulations  to  accentuate  the  steps.  The  first  rise  in  the  profile  ends 
with  the  HEL  associated  with  the  elastic-plastic  transition.  The  second  rise  ends  with  the  volume 
change  from  the  polyamorphic  transformation,  and  the  third  rise  terminates  with  the  maximum 
particle  velocity.  From  this  comparison,  it  is  evident  that  the  elastic-plastic  and  polyamorphic 
transformations  must  occur  concurrently,  or  an  extra  step  will  appear  in  the  velocity  profile. 
Physically,  this  implies  that  the  atoms  move  to  accommodate  shear  strain  as  they  are  rearranging 
to  accommodate  the  volume  change. 
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Figure  62.  Velocity  profiles  from:  (a)  gas  gun  experiments  on  borosilicate  glass 
from  Alexander  et  al.  (23)  and  (b)  simulations  using  independent 
functions  for  the  elastic-plastic  transition  and  the  polyamorphic 
transformation. 


The  modeling  goal  is  to  couple  inelastic  deviatoric  deformation  and  volume  change.  There  is  no 
experimental  data  to  guide  the  functional  form  for  the  coupling,  so  a  simple  quadratic  flow 
potential  model  is  assumed: 


0=  0  =\j aa2e  +  p2  -  p( A). 


(27) 


Here,  ae  is  the  effective  stress,  p  is  the  pressure,  and  p  is  a  function  of  the  volume  fraction  of  the 
low  pressure  polyamorph,  A.  ct  is  a  parameter  regulating  the  relative  influence  of  the  deviatoric 
stress  compared  to  the  pressure. 

Following  concepts  used  in  metal  plasticity,  the  direction  of  inelastic  flow  is  assumed  to  be 
normal  to  the  flow  potential: 


i  _  *30 
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(28) 


where  d,ne/as  is  the  inelastic  part  of  the  rate  of  deformation  tensor,  d,  which  is  decomposed, 
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additively,  into  elastic  and  inelastic  parts  d  =  d eias  +  dme/().,.  The  elastic  part  is  also  a  function  of 
the  stress  tensor  through  the  elastic  moduli.  Together  these  provide  a  coupled  set  of  equations  for 
A  in  equation  28,  which  is  determined  by  an  iterative  procedure. 

The  volume  fraction  of  the  low  density  polyamorph  is  assumed  to  evolve  through  a  kinetic 
relation 


A  =  /3(X 


T  Tre.f  \ 

Tref 


(29) 


where  f3  and  v  are  parameters  and  Aoo  represents  the  equilibrium  fraction  of  the  low  density 
polyamorph.  The  equilibrium  phase  fraction 
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is  a  function  of  the  applied  loading,  p,  the  threshold  pressure  at  which  the  transformation  begins, 
pt ,  and  the  pressure  at  which  the  glass  is  fully  transformed  to  the  high  density  polyamorph,  pf. 

These  equations  were  implemented  into  a  finite  element  code  and  the  gas  gun  simulations  run. 
Figure  63a  shows  the  results  with  the  kinetic  parameter  set  high  (/3  =  100)  to  give  the  equilibrium 
response.  From  these  calculations  it  is  evident  that  a  two-wave  structure  is  produced,  indicating 
that  the  inelastic  deformation  and  polyamorph  transformation  are  occurring  simultaneously. 
Figure  63b  shows  the  results  with  a  kinetic  parameter  (/ 3  =  3)  set  to  reproduce  the  basic  features 
observed  in  the  experimental  data  plotted  in  figure  62a. 

The  simulation  results  demonstrate  some  general  model  features  necessary  to  reproduce 
experimental  observations  from  high-pressure  gas  gun  experiments.  The  model  must  couple  the 
deviatoric  stress  and  the  pressure  for  the  polyamorphic  transformation,  and  a  kinetic  model  is 
necessary  to  introduce  time-dependent  rises  in  the  response.  Other  features  may  also  be 
necessary,  but  these  two  have  a  major  impact  and  provide  a  well-defined  starting  point  for  model 
improvement.  The  particular  functional  relations  assumed  for  this  model  are  conjecture  guided 
by  experience  in  modeling  metals.  Hence,  appropriate  data  are  needed  to  guide  model 
development.  While  it  will  be  possible  to  tune  a  kinetic  parameter  based  on  data  from  current 
experimental  methods,  existing  experimental  techniques  will  provide  little  information  on  the 
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Figure  63.  Gas  gun  simulation  results  for  a  model  with  concurrent  elastic-plastic 
transition  and  polyamorphic  transformation  and  with  (a)  a  high  kinetic 
parameter  for  rapid  transformation  and  (b)  the  kinetic  parameter 
adjusted  to  resemble  experiments. 


functional  form  for  the  kinetic  relation  or  the  flow  function  dependence  on  deviatoric  stress  and 
pressure.  The  alternative,  and  the  approach  pursued  in  the  program,  is  to  determine  the 
functional  relationships  and  parameters  through  multiscale  modeling  approaches. 

9.1  Summary  of  Continuum  Equation-of-state  Model  Development 

The  EOS  development  included  the  non-monotonic  change  in  bulk  modulus  with  pressure  and  a 
polyamorphic  transformation  from  a  lower  density  amorphous  state  to  a  higher  density 
configuration.  The  EOS  form  near  the  reference  density  was  guided  by  previously  published 
pressure-density  data.  Figure  64  depicts  the  model  dependence  of  the  pressure  and  bulk  modulus 
on  density  along  with  the  data  of  Zha  et  al.  (24).  The  agreement  is  reasonable  up  to  10  GPa, 
which  is  near  the  pressure  where  the  polyamorphic  transformation  begins.  The  EOS  for  the  high 
density  polyamorph  was  assumed  to  follow  a  Murnaghan  form  and  the  model  constants  were 
determined  by  fitting  Hugoniot  data  from  Sugiura  et  al.  (67)  and  Alexander  (23).  This  fit  is 
complicated  by  a  lack  of  data.  In  particular,  the  reference  density  for  the  high  density  phase  is 
not  easily  obtained,  since  it  must  be  measured  after  release  from  high  pressures. 

The  remaining  component  of  the  EOS  model  is  the  temperature-dependent  kinetic  relation 
connecting  the  low  and  high  density  polyamorphs.  Here,  it  is  assumed  that  the  transformation 
rate  is  proportional  to  the  cube  of  the  difference  between  the  current  phase  fraction  and  the 
equilibrium  phase  fraction.  The  model  is  assessed  by  comparing  interface  velocity  predictions 
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Figure  64.  Model  fit  of  pressure  and  bulk  modulus  with  data  from  Zha  et  al.  (24). 


from  gas  gun  simulations  with  available  data  in  figure  62a.  The  curvature  of  the  velocity  data  is 
due  entirely  to  the  kinetic  relation,  and  it  is  capturing  much  of  the  inflection  and  curvature.  The 
timing  of  the  rises  in  the  curves  is  off,  which  suggests  that  the  underlying  response  may  be  too 
stiff  in  the  elastic  range  and  that  the  volume  change  associated  with  the  transformation  may  be  too 
large.  These  can  be  adjusted  when  more  data  are  available. 


10.  A  Perfectly  Matched  Layer  for  Peridynamics  in  Two  Dimensions 


In  this  section,  we  develop  a  peridynamic  method  for  modeling  nanoindentation  experiments  in 
fused  silica  and  other  chemically  substituted  glasses.  Originally  introduced  in  Silling  (68), 
peridynamics  is  a  non-local  formulation  of  elastodynamics,  which  can  more  easily  incorporate 
discontinuities  such  as  cracks  and  damage.  Derivatives  of  field  variables  in  the  classical 
continuum  model  are  replaced  by  integrals  over  a  small  neighborhood  of  microelastic  kernels, 
which  replaces  the  standard  constitutive  relation.  In  its  discretized  form,  an  elastic  solid  is 
treated  as  a  collection  of  particles  or  nodes,  each  connected  to  its  neighbors  by  breakable  bonds. 
Bond  breakage  can  be  defined  to  occur  when  a  bond  is  stretched  past  some  predetermined  limit. 
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The  end  result  is  a  method  capable  of  predicting  crack  growth  in  brittle  elastic  materials  (69-74). 
The  original  formulation  was  limited  to  materials  with  a  fixed  Poisson’s  ratio  of  0.25;  however, 
state-based  peridynamics  was  introduced  allowing  for  more  flexible  constitutive  relations  (75). 

As  becomes  clear  later,  state-based  peridynamics  allows  for  an  auxiliary  field  formulation,  which 
is  necessary  to  implement  a  perfectly  matched  layer  (PML).  While  most  peridynamics  work  has 
focused  on  simulating  problems  with  free  or  fixed  boundary  conditions,  there  are  applications  in 
which  the  simulation  of  an  infinite  medium  may  be  useful,  such  as  crack  propagation  in  a 
halfspace  or  nanoindentation  problems.  Absorbing  boundary  conditions  are  one  way  of 
simulating  an  infinite  medium  as  any  impinging  waves  are  suppressed  so  they  do  not  reflect  back 
into  the  simulation.  A  PML  is  such  an  absorbing  boundary,  and  was  originally  introduced  for 
electromagnetic  simulations  (76,  77).  PMLs  differ  from  traditional  absorbing  boundary 
conditions  in  that  they  are  an  absorbing  layer  with  a  finite  width,  placed  between  the  computation 
region  of  interest  and  the  truncation  of  the  grid  or  mesh.  They  can  also  be  thought  of  as  an 
anisotropic  absorbing  material,  which  is  why  the  flexibility  of  a  state-based  peridynamics  is 
necessary.  A  PML  was  applied  to  one-dimensional  (1-D)  peridynamics  (78),  which  used  the 
results  of  Du  et  al.  (79)  to  formulate  an  auxiliary  field  equation.  This  approach  required  a  matrix 
representation  of  the  auxiliary  field,  which  may  be  memory  prohibitive  in  higher  dimensions. 

10.1  Two-dimensional  (2-D),  State-based  Peridynamics 

The  continuum  equation  of  motion  in  an  elastic  solid  can  be  stated  as 

d2  _ 

p— u  =  V  •  er  +  b  =  V  •  (c  :  e)  +  b,  (31) 

where  p(x)  is  the  density;  u(x,  t)  is  the  displacement;  <r(x,  t)  is  the  stress  tensor;  e(x,  t)  is  the 
strain  tensor;  c  is  the  stiffness  matrix  for  plane  strain,  defined  by  Young’s  modulus  E,  Poisson’s 
ratio  v,  or  Lame  parameters  A  and  //;  and  b(x)  is  a  body  force  (80).  (Throughout,  boldface  type 
denotes  a  vector  and  a  boldface  variable  with  an  overbar  denotes  a  tensor.)  Equation  31  is  a  local 
formulation  because  the  divergence  of  the  stress  (and  gradient  of  the  displacement  implied  in  its 
definition)  represents  a  local  operation  on  a  variable.  In  other  words,  the  action  of  V  •  rr  only 
depends  on  er  at  a  single  spatial  point.  In  problems  involving  discontinuities,  such  as  cracks,  the 
divergence  at  such  discontinuities  is  not  well  defined,  leading  to  numerical  implementation 
problems.  Peridynamics  proposes  replacing  V  •  rr  with  a  nonlocal  operation  that  nonetheless 
also  represents  a  force.  Here,  we  use  the  state-based  peridynamics  (75),  rather  than  the  original 
bond-based  version  (68).  A  PML  application  requires  an  auxiliary  field  formulation,  as  it  is 
essentially  an  anisotropic  absorbing  material,  if  a  non-physical  one.  Consequently,  a  state-based 
peridynamic  formulation  is  necessary  to  implement  the  required  constitutive  relations  in  the 
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absorber.  The  state -based  peridynamics  uses  a  family  of  bonds  to  determine  a  given  force  rather 
than  a  single  bond  independently.  This  more  general  approach  allows  for  inelastic  behavior  and 
more  general  elastic  behavior,  and  is  governed  by 

d2  f  —  — 

Pgpu=  (T[x,  t] (x'  —  x)  —  T[x',  t] (x  —  x'))  dV1  +  b,  (32) 

where  Hx  is  the  horizon  region,  defined  as  a  circle  centered  at  x  with  radius  5,  and  T[x,  t]  (x'  —  x) 
is  a  peridynamic  vector  state,  with  parameters  in  the  square  brackets  indicating  variables  that  act 
as  arguments  to  any  functions  that  define  the  vector  state  and  variables  in  the  angle  brackets 
acting  as  arguments  to  the  vector  state  itself.  In  the  state -based  formulation,  the  deformation 
gradient,  given  by 


F  =  I  +  uV, 


(33) 


can  be  approximated  as  a  vector  state  as 


F[x,  t]  = 


Cm(Y[x,t}(0®Z)dVxl 


UHx 


IT1, 


(34) 


where  C(|£|)  =  exp(—  |£|2  / 6 2)  is  a  shape  function,  taken  as  a  Gaussian  distribution  here  and 
with  the  horizon  H.x  extended  so  that  the  shape  function  decays  to  an  arbitrary,  small  value,  taken 
here  as  10-6,  K  is  a  shape  tensor  given  by 


K[x,  t]  —  /  £(!£!)  (£<8>|)dK,,K 
Jux 

and  Y  is  a  deformation  vector  state  given  by 


-l 


iLinv  uim 
^ xx  ^ xy 

b.  inv  h,inv 
^ yx  ^ yy 


(35) 


Y[x,f](£)  =  rj  +  £, 


(36) 


with  r]  =  u[x',  t]  —  u[x,  t]  and  $,  =  x'  —  x  (81). 

The  deformation  gradient  can  now  be  substituted  into  Hooke’s  law  and  strain-displacement 
relations,  giving  a  stress  term  W  in  terms  of  u  in  plane  strain 

d2  _ 

p—u  =  V  •  cr  =  V  •  (c  :  e)  , 


(37) 
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where 


e[x,  t]  =  ^  (Vu  +  uV)  «  i  (F[x,  t]  +  F[x,  t]T  -  21)  . 
Ultimately,  the  peridynamic  vector  state  T  for  plane  strain  elasticity  is  given  by 

T[x,m=C(\ZM*,t]K1£. 


(38) 


(39) 


10.2  Auxiliary  Field  Formulation  and  PML  Application 


The  first  step  in  formulating  a  PML  is  to  construct  an  analytic  continuation  to  the  complex  plane, 
such  as  x  =  x  +  ig(x),  where  g(x)  is  a  given  function  describing  the  deformation  (82).  This 
mapping  has  the  effect  of  transforming  traveling  waves  of  the  form  elkx,  where  k  =  u/c  is  the 
wave  number,  into  evanescent  waves  of  the  form  elkxe~k9(x)  ^  thus  attenuating  such  waves  in  the 
PML  region.  Applying  a  PML  involves  substituting  for  spatial  derivatives  using 


d  Id 

dx  i  +  dx 

UJ 


(40) 


The  function  cp(x)  defines  the  extent  of  the  PML  region,  i.e.,  when  cf>(x)  =  0  the  original  wave 
equation  is  obtained,  and  when  <f>(x)  >  0,  traveling  waves  decay  exponentially.  Typically,  <p(x) 
transitions  from  0  to  a  constant  value  using  a  smooth  function  to  prevent  numerical  reflections  at 
the  boundary  between  the  absorbing  and  computational  regions.  Before  applying  a  PML  directly 
to  the  peridynamic  equation,  equation  3 1  is  treated  so  that  the  PML  application  to  peridynamics  is 
clear.  It  is  convenient  to  convert  equation  3 1  to  the  Laplace  domain,  assuming  e~st  time 
dependence,  and  express  the  wave  equation  as  two  coupled  first-order  partial  differential 
equations,  the  first  in  u  and  the  second  in  =  W 


psu  —  V  •  if} 
sip  —  c  :  e, 


(41) 


where  the  Laplace  transform  of  a  variable  is  indicated  by  £{/}  =  /.  Expanding  equation  41  into 
components  gives  five  coupled  equations 
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psux 

pSUy 

si>x 

S'ljjy 

si>T 


d  ~  d  ~ 

+  7 T~Wr 

ox  ay 

d  ~  d  ~ 

+ 

c\  \  d  ~  .  d  „ 

(A  +  2 ^  dxUx  +  XdyUy 


d 


d 


dx~x  '  v"  '  d y~y 


^~nZUx  +  (-^  +  2 p)  -^—U 

P 


d_ 

dy 


d_ 

faUy 


(42) 


We  can  expand  the  state-based  formulation  into  components  and  match  terms  to  equation  42. 
Following  this  approach  yields  a  viable  method  for  performing  PML  substitutions.  First,  the 
state-based  peridynamic  equations  listed  above  in  equations  32-39  can  be  written  explicitly  as 


psux[x,  s}=  C(|£|)  ( ipx[x,  s\k xx  +  ij>T[x,  s\k™ )  £x  +  ( $x[x',  sjk ™  +  Vv[x',  s]k™  )  £ 


C(l«l)  #1”  +  &[*,  #”v)  „  +  («*',  s]k‘"'  +  P[x’,  »]*£')  {, 


iHx 


dVx> 
dVx>, 
dVx> 

C'(l^l)  \(^r[x,s}k™  +  ^y[x,s]k™)£y+  UT[x!,s]V™' +  i>y[x!,s]V™')£y\  dVx, , 


psuy[x,  s}=  C(|£|)  |  ( i/>T[x,  s]k™  +  £>y [x,  s]k% )  £x  +  ( i/>T[x!,  s]k™  +  ipy[x',  s]k™  )  £x 

Jnx 


si>x[x,s]  =  (A  +  2  p) 


Uhx 


C'dCl)  Yx[x,  s]£xk ™  +  Yx[x,  s]£yk™  dVx,  -  1 


+  A 

s$y\x,  s]  =  A 


Uh* 


<7(1*1)  Yy[x,  s}£xk™  +  Yv[x,  s}£yk™  dVx,  -  1 


(43) 


(A  +  2/i) 


<7(1*1)  [Yx[x,  s]£xk™  +  Yx[x,  s]£yk dVx,  -  1 

C-(ICI)  fe[*.  0%kZ  +  ?„[*,  «]{/•"■)  iV*  -  1 


>HX 


s$T[x,  s}=  p  C'(lll)  [Yx[x,  s\£xk™  +  Yx[x,  s}£yk™)  dVx> 


+  P  /  m\)(Yylx,s]£xk™  +  Yy[x,s]£yk™)dVx., 
Jnx 


etc.  Though  no  derivatives  appear  in  equations  43,  the  correspondence  of  each  term  to  those  in 
equation  42  is  apparent  and  the  PML  substitutions  can  be  made.  For  example,  the  first  equation 
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in  43  can  be  rewritten  as 


P  (s  +  (fix)  (s  +  (fiy)  Ux  =  (s  +  (fiy)  /  C(|£|) 

Jh x 

+  (s  +  (fiy )  [  C'(IH) 

JHx 

+  (s  +  (fix)  [  C(|£|) 

Jhx 

+  (s  +  (fix)  [  C(|*|) 

JUy. 

with  the  remaining  equations  following  similarly, 
to  the  time  domain  and  implementing  a  forward  Euler  discretization  scheme  in  time,  and  the 
standard  one-point  integration  method  (74)  in  space. 

10.3  Results 

The  previous  method  was  numerically  implemented  using  a  forward  Euler  method  for  the 
temporal  discretization,  and  standard  one-point  integration  with  point-matching  in  space.  The 
PML  was  tested  on  three  types  of  problems:  a  wave  propagation  problem  to  demonstrate  the 
effectiveness  of  the  PML,  two  crack  propagation  problems  (one  state -based  and  one  bond-based), 
and  an  axisymmetric  indentation  problem. 

10.3.1  Wave  Propagation 

The  PML  was  first  tested  on  a  wave  propagation  problem  with  PML  boundary  layers  and  a 
Gaussian  distribution  as  an  initial  condition.  Specifically,  the  x-di  reeled  displacement  was  set  to 

ux(x,t  =  0)  =  200lx— Pmidl2  ^  (45) 

where  pmid  is  the  mid-point  of  the  region,  which  in  this  example  was  defined  as  0  <  x,  y  <  1  and 
discretized  with  Ax  =  Ay  =  0.01  m.  Young’s  modulus  for  the  region  was  set  to  1  Pa,  Poisson’s 
ratio  was  1/4,  and  the  density  was  1  kg/m3.  The  PML  region  was  defined  as  the  0.3-m  border 
around  the  1  m-by-1  m  region  and  used  a  Gaussian  ramp  with  a  width  of  0.2  m,  finally  reaching  a 
maximum  of  50  s_1  for  the  remaining  0.1  m.  Lor  the  Gaussian  kernel,  a  horizon  size  of 
5  =  l.lAx  was  used,  and  for  the  Heaviside  kernel,  a  horizon  of  5  =  3.1  Ax  was  used. 

The  simulation  was  run  with  both  the  Heaviside  and  Gaussian  kernels,  with  the  total  strain  energy 
shown  in  figure  65.  The  Gaussian  kernel  (dotted  line)  shows  the  largest  drop  in  energy,  reaching 
a  minimum  of  5.6  x  10-7,  and  the  Heaviside  kernel  (dashed  line)  decreases  to  1.  x  10-4. 


^*[X,  SWxx  +  i’xyi*,  S}k™  )  CxdVy 


^[x',  s]k%'  +  $xvW,  j  txdV, 


(44) 


ifix[x,  s]k™  +  ifiT[x,  s] k'yyj  £ydV* 


The  final  step  involves  simply  converting  back 
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Figure  65.  Total  strain  energy  in  a  simulation  terminated  by  a  PML. 


A  bounded  simulation  is  shown  for  reference  as  the  solid  line,  which  used  a  fixed  displacement 
boundary  condition  and  the  Gaussian  kernel.  Figure  66  shows  a  waterfall  plot  of  the  x-dircctcd 
displacement  along  the  y  =  0.5  line  for  the  Gaussian  kernel  with  the  PML  function  (j>x  shown  in 
gray  on  the  far  end  of  the  plot  (corresponding  to  t  =  1  s).  Figure  67  shows  the  absolute  value  of 
the  x-dircctcd  displacement  at  the  edge  of  the  PML  region,  in  simulations  terminated  by  a  PML 
and  with  a  fixed  boundary  condition.  The  wave  is  absorbed  at  the  boundary  with  minimal 
reflections:  as  can  be  seen,  the  plots  align  for  a  time,  and  where  they  deviate  (indicating  a 
reflection  from  the  hard  boundary),  the  PML  simulation  remains  in  decay. 


For  verification,  the  method  was  compared  with  an  exact  analytical  solution.  Consider  a 
cylindrically  symmetric  wave  propagating  in  an  infinite  elastic  medium  with  the  same  constitutive 
parameters  as  the  above  example,  and  with  an  initial  condition  given  by 


“»(r)  = 6  © 


•f  \  2 

l+(- 

a 


(46) 
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(b  =50  s'1 
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Figure  66.  x-directed  displacement  at  y  =  0.5  m,  terminated  by  a  PML. 


Figure  67.  x-directed  displacement  at  x  =  0.3,  y  =  0.5  m.  The  solid  line  shows 
results  terminated  by  a  PML,  and  the  dashed  line  used  a  fixed 
boundary  condition. 
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where  here  we  take  6  =  1  and  a  =  0.1.  The  exact  solution  is  given  by  Eringen  and  Suhubi  (83) 


br 


u(r, *)  =  7k  nry/R2  +  a  (2«  -  R2) , 


O  —  1  “h 


V2aR 6 
r2  —  c2t2 


R2  =  Ja2  + 


4c2t2 


(47) 


where  c  is  the  longitudinal  wave  speed.  This  problem  was  simulated  in  a  2-D  region,  2  m-by-2  m 
and  Ax  =  Ay  =  0.01  m,  terminated  by  a  PML  with  the  same  dimensions  and  magnitude  as  the 
above  problem.  The  Gaussian  kernel  was  used  with  a  horizon  size  of  5  =  0.73  Ax,  with  an  actual 
cutoff  of  0.028  m.  The  results  are  shown  in  figure  68,  with  the  exact  solution  shown  as  the  solid 
line  and  the  peridynamic  solution  shown  as  the  dashed  line.  The  peridynamic  solution  shows 
good  agreement  with  the  exact  solution  and  minimal  reflections  from  the  PML  boundary. 


Figure  68.  x-directed  displacement  at  x  =  0.5,  y  =  0  m.  The  solid  line  shows 
the  exact  solution  and  the  dashed  line  results  terminated  by  a  PML. 


10.3.2  Crack  Propagation 

Crack  propagation  in  a  half-space  can  be  useful  for  modeling  physical  phenomenon  such  as 
indentation  experiments.  As  an  example,  we  model  such  a  problem  as  a  body  force  applied  to  a 
finite  region  with  small  pre-cracks  in  a  region  terminated  on  three  sides  with  PMLs.  One  addition 
to  the  algorithm  for  this  problem  was  a  drag  term,  used  to  reduce  noise.  For  crack  problems  with 
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a  sudden  force  application,  noise  and  oscillations  can  cause  hot  spots  and  undesirable  cracking. 

To  remedy  this,  a  drag  term  can  be  added  to  smooth  oscillations,  by  adjusting  the  nodal  velocity  as 


(48) 


where  D  is  the  drag  coefficient  and  Nb  is  the  remaining  number  of  bonds  in  the  family  of  node  n 
(84). 


An  absorbing  boundary  ensures  that  no  reflections  from  the  boundaries  interfere  with  the  crack 
propagation,  possibly  causing  it  to  deviate.  Figure  69  gives  a  schematic  of  the  problem,  the  extent 
of  the  computation  region  is  designated  by  the  solid  line,  the  PML  ramp  begins  at  the  dashed  line 
and  the  PML  plateaus  at  the  dotted  line.  The  computational  region  was  70  mm  wide  and  35.25 
mm  high,  the  PML  region  began  at  15  mm  from  each  edge  (except  the  top)  and  peaked  at  5  mm  to 
a  value  of  5  x  106  s'1.  The  node  spacing  was  0.496  mm  and  the  time  step  size  was  1  ns.  For 
material  values,  the  density  was  2235  kg/m3,  Young’s  modulus  was  65  GPa,  Poisson’s  ratio  was 
0.2,  and  the  fracture  criteria  used  a  fracture  energy  of  204  J/m2.  The  failure  criteria  used  in  this 
simulation  was  bond-based,  i.e.,  a  bond  failed  if  it  was  stretched  past  a  given  limit,  determined  by 
the  fracture  energy  (71).  The  maximum  relative  bond  stretch  was  then  2.971  x  10  3.  The  load 
was  applied  across  a  10-mm  region,  centered  at  the  top  surface,  with  pre-cracks  on  each  edge  with 
a  length  of  two  nodes  or  0.993  mm.  The  simulation  was  run  for  a  total  of  10  /is,  and  the  cracks 
were  measured  manually  from  the  edge  of  the  pre-cracks  to  the  extent  of  the  damaged  area. 


00 


00 


00 


Figure  69.  Schematic  diagram  of  the  problem  for  crack  propagation  in  a 
half-space. 
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Figure  70  shows  a  result  that  had  an  applied  load  of  250  x  109  N,  yielding  a  2.1 1-mm  crack. 
Figure  71  shows  a  close  up  of  the  damaged  area  from  figure  70.  As  can  be  seen,  the  crack 
extends  three  nodes  down  and  three  nodes  across.  Finally,  the  applied  load  was  varied  between 
140  x  109  N  and  500  x  109  N,  with  the  crack  length  versus  applied  load  shown  in  figure  72  as  the 
dots.  A  curve,  shown  as  the  solid  line  in  figure  72,  was  fit  using  the  form  motivated  by  the 
universal  scaling  law  presented  in  section  6: 

£  =  Aps ,  (49) 

where  i  is  the  crack  length  in  meters  and  p  is  the  applied  load  in  newtons.  A  linear  fit  was 
computed  in  log-log  space  using  least  squares  giving  a  power  law  of  s  =  1.96  and 
A  =  4.3  x  10-8  m/N'\  This  fit  matches  well  with  the  predicted,  dimensional  analysis  law  shown 
in  equation  17.  Though  the  loading  is  different  here,  because  the  problem  is  2-D  the  dimensional 
analysis  still  applies  giving  the  correct  power  law  form.  In  other  words,  for  all  2-D  problems 
with  stable  crack  growth  in  infinite  regions,  we  should  expect  a  load-crack  length  relation  similar 
to  that  of  equation  49  with  a  power  of  s  =  2,  while  in  3-D,  we  expect  a  power  of  s  =  2/3. 


Figure  70.  Damage  map  resulting  from  a  250  x  109  N  applied  load  after  10  //s. 
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Figure  71.  Close  up  of  damage  map  from  figure  70. 


Figure  72.  Crack  length  versus  applied  force  for  indentation  into  an  elastic 

half-space  using  state-based  peridynamics.  Blue  dots  represent  data 
points  from  the  peridynamic  simulation  and  the  solid  line  is  a  curve  fit. 
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10.3.3  Bond-based  Verification  of  a  Center  Crack 


A  center  crack  opening  problem  was  simulated  with  a  bond-based,  non-PML  terminated 
peridynamics  method.  The  problem,  illustrated  schematically  in  figure  41,  involves  applying 
equal  but  oppositely  directed  point  loads  in  an  infinite  region.  The  stable  crack  length  is  related 
to  the  applied  load  as  shown  in  equation  17.  Here,  rather  than  applying  the  load  at  a  single  node 
(which  may  cause  instability  or  large  oscillations),  the  load  was  spread  over  three  nodes  on  each 
side  at  the  center  of  a  square  region.  Bond-based  peridynamics  was  used  with  a  node  spacing  of 
0.5  mm,  a  time  step  size  of  1  ns,  Young’s  modulus  of  65  GPa,  density  2235  kg/m3,  and  horizon 
size  of  1.5  mm  with  a  Heaviside  shape  function.  The  load  was  applied  slowly  with  a  Gaussian 
ramp,  spread  over  500  time  steps.  The  region  was  not  terminated  with  a  PML,  though  the  region 
size  (100  mm-by-100  mm)  was  large  enough  to  avoid  issues  associated  with  crack-boundary 
interactions.  The  criterion  for  bond  damage  was  a  critical  relative  stretch  of  0.00185.  A  small 
pre-crack  was  added  between  the  two  point  loads  and  the  bonds  surrounding  the  loads  were  set  to 
a  no-fail  condition.  This  setup  was  run  for  varying  loads,  between  4  x  1011  and  9  x  1011  N  and 
the  crack  half  length  was  measured  between  the  center  of  the  region  and  the  crack  tip,  with  a 
crack  tip  being  a  region  with  a  broken  bond  ratio  of  more  than  30%.  The  results  are  shown  in 
figure  73  on  a  log-log  scale.  As  in  the  above  example  (see  equation  49),  a  linear  curve  fit  in 
log-log  space  (shown  in  figure  73  as  the  solid  red  line)  was  computed  using  least  squares  to 
determine  the  exponent  of  the  power  law,  and  found  to  be  1 .94,  very  near  2,  the  value  determined 
by  dimensional  analysis  in  equation  41. 


Figure  73.  Crack  half  length  versus  load  for  a  center  crack  problem  using 
bond-based  peridynamics. 
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10.3.4  Verification  of  an  Axisymmetric  Indentation  Problem 


Finally,  the  state-based,  PML  formulation  of  peridynamics  was  implemented  in  an  axisymmetric 
setting,  the  details  of  which  will  be  published  in  the  future.  Sneddon  (25)  derives  an  exact 
analytical  solution  for  a  flat,  rigid,  cylindrical  indentor  being  pressed  into  an  infinite  elastic 
halfspace.  Given  the  elastic  properties  of  the  medium,  the  depth  of  indentation,  and  the  indentor 
radius,  the  displacement  and  stress  can  be  computed  for  the  entire  region.  To  simulate  this 
problem,  an  infinite  region  was  represented  using  a  PML  in  a  finite  region  with  a  radius  of  50  mm 
and  a  depth  of  50  mm.  The  flat  extent  of  the  PML  was  5  mm  and  the  ramp  region  was  15  mm. 
The  remaining  parameters  were  identical  to  those  in  the  above  examples,  though  no  bond 
breaking  was  used  in  this  example.  The  indentor  had  a  radius  of  5  mm,  and  was  pressed  to  a 
depth  of  10  //m.  Sneddon’s  exact  solution  (25)  for  the  ^-directed  displacement  at  the  surface  of 
the  halfspace  is  given  by 


uz  ( r ,  z  =  0) 


—  sin  1  (-)  .  r  >  a 

7r  V  r  )  1 

e,  r  <  a 


(50) 


where  e  is  the  indentation  depth  and  a  is  the  radius  of  the  indentor.  The  results  of  the  simulation 
are  plotted  in  figure  74  with  the  peridynamic  solution  shown  in  blue  and  the  exact  solution  in  red. 
In  addition,  the  magnitude  of  the  PML  is  indicated  as  the  graded  gray  region  on  the  right  of  the 
figure,  labeled  as  “PML.”  As  can  be  seen,  the  peridynamic  solution  is  forced  to  zero  in  this 
region,  and  is  not  expected  to  match  the  exact  solution.  The  main  discrepancy  is  due  to  the 
discontinuous  load  applied,  which  leads  to  ringing  in  the  peridynamic  solution.  Remedies  for 
this  issue  will  be  researched  in  the  future. 
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Figure  74.  Peridynamic  solution  (blue)  of  an  indentation  problem  compared  with 
Sneddon’s  exact  solution  (25)  (equation  50)  (red). 


11.  Glass  Technology  Short  Course 


A  one-day  short  course  on  the  “Fundamentals  of  Glass  Science,”  was  taught  by  Professor  Arun  K. 
Varshneya,  Alfred  University,  at  ARL  on  October  29,  2010.  Course  attendees  received  copies  of 
Varshneya’s  Fundamentals  of  Inorganic  Glasses  (26).  The  course  covered  topics  on  (1)  basic 
compositions  and  families  of  commercial  glasses  that  include  vitreous  silicates,  SL  silicates, 
borosilicates,  lead  silicates,  and  aluminosilicates;  (2)  fundamentals  of  the  glassy  state;  (3)  phase 
separation  and  liquid-liquid  immiscibility;  (4)  glass  structures  (oxide  and  non-oxide  glasses);  (5) 
review  of  some  glass  properties  (e.g.,  elastic  properties,  hardness,  viscosity,  thermal  expansion, 
durability);  and  (6)  annealing  and  strengthening. 

The  course  was  well  received,  and  a  question-and-answer  period  ensued  with  questions  like, 
“What  is  the  difference  between  a  glass  and  an  amorphous  material?”  Answer:  “Amorphous 
materials  make  a  continuous  or  discontinuous  volume  transition  to  a  crystal  or  vapor  on  heating, 
whereas,  glasses  continuously  change  to  a  liquid  on  heating.”  See  also,  p.  17  and  18  in 
Varshneya  (26),  and  figure  75. 
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Figure  75.  The  volume-temperature  diagram  for  a  glass-forming  liquid,  after 
Varshneya  (26). 


12.  A  Short-term  Conceptual  Project 


A  short-term  conceptual  project  to  determine  an  effective  experimental  and  theoretical  approach 
to  model  and  characterize  the  role  of  glassy  materials  in  resisting  ballistic  impact  was  conducted 
by  Richard  Lehman,  Professor  and  Chair  of  the  Department  of  Materials  Science  and 
Engineering,  Rutgers  University,  Piscataway,  NJ.  A  short  synopsis  of  the  report  was  briefed  at 
ARL  on  November  17,  2011,  and  appears  below. 

A  number  of  features  of  the  glassy  state  are  thought  to  participate  in  the  ballistic  response  of  glass: 

•  Structural  relaxation/viscoelasticity/strain  energy  equivalence  theory  (SEET) 

•  Short-  and  longer-range  atomic  structural  characteristics 

•  Free  volume  (compaction,  bulking) 

•  Cation  and  anion  coordination 

•  Bond  energies 
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•  Other  nanoscale  order  characteristics  that  are  difficult  to  determine  experimentally  for  silicate 
glasses. 


A  number  of  expert  personnel  were  contacted  with  the  following  panel  discussion  outcomes: 


•  Structural  relaxation:  The  ability  of  glass  to  convert  to  a  liquid  with  little  structural 
adjustment  compared  to  crystalline  materials  may  enable  glass  to  respond  more  favorably  to 
shaped-charge  assault  within  the  microsecond  time  periods  of  the  impact.  Viscoelastic  issues 
will  be  addressed. 

•  Modeling:  Structural  models  based  on  statistical  thermodynamics,  ring  structures,  molecular 
dynamics,  and,  anisotropic  finite  element  methods  (FEMs)  will  be  considered.  Experimental 
pair  distribution  functions  based  on  x-ray  data  will  be  discussed. 

•  Molecular  defects:  Glass  is  rich  in  molecular  defects  such  as  oxygen  hole  centers,  E, 
bridging  oxygen,  non-bridging  oxygen,  and  various  ring  defect  structure.  These  defects  may 
play  an  important  role  in  the  shaped-charge  behavior. 

•  Heterogeneous  structure:  Most  glasses  are  not  the  uniform  isotropic  material  at  the 
mesoscale  that  many  scientists  assume.  Phase  separation,  cation  clustering,  and  small-scale 
density  fluctuations  (as  evidenced  by  various  types  of  optical  scattering)  may  be  important  in 
allowing  glass  to  accommodate  high  strain  rates. 

•  Free- volume:  Glasses  contain  a  large  amount  of  free  volume  compared  to  their  crystalline 
counterparts. 


Summary  and  conclusions  of  the  study  include  the  following: 


•  MDs  can  be  expanded  into  the  mesoscale  region: 

-  Large  atom  arrays  (>109)  and  non-cubic  shapes,  e.g.,  high  aspect  ratio  cylinders. 

-  Use  surrogates  to  boost  scale. 

-  Time  scale  (microseconds)  cannot  be  collapsed  and  is  a  major  limit  on  modeling  time. 

•  FEMs  can  be  down-scaled  to  the  mesoscale  size  region: 

-  Requires  specialized  nonlinear  methods. 

-  Structural  elements  modeled  as  nonlinear  elements. 

-  Xi  Chen  at  Columbia  Univ.  has  relevant  experience. 
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•  Medium-ranged  structural  data  obtained  with  scanning  tunneling  electron  microscopy 
(STEM): 

-  Initial  material  characterization. 

-  Input  to  modeling  effort. 

-  Aid  in  interpretation  of  stress  wave  work. 

•  Stress  wave  characterization: 

-  Coherent  acoustic  phonon  spectroscopy  is  a  promising  dynamic  approach. 

-  The  time  scale  is  very  short. 

-  Experiments  can  be  adjusted  to  span  a  range  of  experimental  time  periods. 


13.  Effect  of  Glass  Composition  on  the  Performance  of  Glass:  Understand 
the  Role  of  Boron  Concentration  on  the  Dynamic  Properties  in  the 
Borosilicate  System 


Glass  is  currently  used  for  transparent  armor  systems  and  of  interest  for  opaque  armor 
applications.  The  glasses  that  have  been  studied  in  previous  investigations  have  been  commercial 
glasses  such  as  SL  silica  glass,  borosilicate  glass,  and  fused  silica.  Attempts  to  procure  new  glass 
composition  have  been  limited  to  bulk  glasses  for  compositions  that  have  a  commercial  market. 
This  is  primarily  due  to  the  large  capacity  nature  of  glass  production.  Typical  melts  range  from 
small  volume  of  less  than  1  ton  a  day  to  very  large  at  greater  than  100  tons  a  day.  These  are  all 
greater  than  the  needs  required  for  small-scale  evaluations  for  ballistic  and  dynamic  properties. 
Thus,  a  facility  at  ARL-WMRD  to  produce  small  boules  of  glass  to  conduct 
composition-structure-property  relationship  studies  was  installed  in  fiscal  year  2013. 

Figure  76  is  a  photograph  of  a  bottom  loading  furnace  located  at  ARL-WMRD.  It  has  been  used 
for  producing  SL  silica  glass  boules  as  shown  on  figure  77.  These  glasses  were  produced  to 
understand  the  temperature  profile  and  capability  of  the  furnace.  The  fiscal  year  2013  effort  will 
investigate  the  role  of  boron  concentration  on  the  properties  of  glass.  The  approach  is  using  glass 
frits  to  cast  round  disks  of  borosilicate  glasses.  The  target  compositions  are  the  commercial 
Borofloat  composition,  the  Borofloat  composition  with  an  additional  5  wt.%  and  the  Borofloat 
composition  with  a  reduced  5  wt.%  of  boron.  Chemical  analysis,  quasi-static  property 
characterization,  dynamic  properties  characterization,  and  ballistic  performance  will  be  measured 
and  reported. 
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Figure  76.  Bottom  loading  furnace  installed  at  ARL. 


Figure  77.  Typical  glass  boule  produced  in  a  silica  crucible. 
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14.  Conclusions 


This  third-year  final  report  on  multiscale  modeling  of  noncrystalline  ceramics  (glass)  has  focused 
on  establishing  the  framework  for  development  of  a  multiscale  computational  methodology  for 
optimizing  or  enhancing  the  performance  of  silicate-based  glass  materials  not  yet  synthesized.  A 
more  immediate  research  objective  is  to  understand  why  certain  chemically  substituted 
silica-based  materials  exhibit  enhanced  performance  in  the  defeat  of  SCJs  and  other  ballistic 
threats.  Conclusions  consistent  with  the  milestones  shown  in  the  five-year  roadmap  shown  in 
figure  3  are  as  follows: 

1.  Fused  silica  and  various  chemically  substituted  silicate-based  specimens  were  delivered  to 
ARL  under  the  auspices  of  an  ongoing  cooperative  agreement;  these  specimens  have  been 
ballistically  tested,  and  serve  to  validate  future  computational  models  of  fused  silica. 

2.  The  structural  characteristics  of  silicate-based  glasses  include  SRO  consisting  of  atom  to  atom 
bond  lengths  that  are  less  than  0.5  nm  and  bond  angles  characterized  by  RDFs  and  SAXS 
measurements.  IRO  in  silica-based  glasses  consists  of  the  polymerization  of  the  silica  atomic 
tetrahedra  (one  Si  atom  surrounded  by  four  0  atoms  into  ring  structures  of  joined  tetrahedra 
of  a  variety  of  sizes  ranging  from  four  to  eight  groups  of  tetrahedra).  Chemical  substitution 
of  Na,  K,  Mg ,  Co ,  B.  and  other  atoms  can  have  a  profound  effect  on  IRO  and  ballistic 
performance.  Initial  free  volume  and  its  variation  with  pressure  control  densification 
behavior,  but  this  is  difficult  to  measure  or  quantify,  especially  at  elevated  pressure.  A  glass 
brittleness  parameter  discussed  by  Ito  (4)  appears  dependent  on  the  deformation  and  fracture 
behavior,  which  depends  on  flow  and  densification  before  crack  initiation  and  on  the  bond 
strength  of  the  network  and  seems  to  decrease  with  a  decrease  in  density. 

3.  A  series  of  experiments  on  Borofloat,  Starphire,  and  fused  silica  were  conducted  at  the 
Emst-Mach  Institute  in  Germany,  which  showed  that  fracture  and  fragmentation  (figure  21)  of 
these  glasses  have  profoundly  different  macroscopic  response  (figure  12)  and  fracture  kinetics 
(table  8)  useful  for  computational  model  validation.  Solid  cylinder  impact  onto  fused  silica 
targets  results  in  a  greater  mass  of  finer  fragments  than  the  AP  round,  whereas  the  AP  round 
results  in  a  larger  mass  of  greater  than  2-mm-size  fragments  (figure  21). 

4.  A  series  of  nanoindentation  experiments  (indenter  radius  =  3  //m)  into  fused  silica  were 
conducted  at  ARL,  which  showed  that  both  hardness  and  modulus  decreased  with  increasing 
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indentation  depth  (figure  30).  Also,  both  radial  and  cone  cracks  were  observed  in  indents, 
which  exceeded  1500  nm  (figures  32  and  33).  The  nanoindentation  experiments  will  be  used 
to  validate  MD  and  peridynamics  models  of  fused  silica  and  other  chemically  substituted 
glasses. 

5.  High-pressure  DAC  experiments  were  conducted  at  ARL  to  pressures  of  ~33  GPa  and 
provide  a  means  to  measure  density  changes  with  pressure  in  glasses.  Our  DAC  results  on 
fused  silica  compared  well  with  prior  published  results  (72)  and  provide  a  means  for 
validation  of  AIMD-based  predictions  of  Hugoniots  for  more  complex  chemically  substituted 
glasses.  Preliminary  measurements  of  neutron  diffraction  spectra  of  fused  silica,  Borofloat, 
and  SL  glasses  provide  an  independent  means  to  assess  changes  in  glass  network  structure 
with  pressure.  Preliminary  analysis  of  Raman  spectroscopic  measurements  on  several  glasses 
with  various  (known)  concentrations  of  secondary  molecules  reveal  changes  in  the  300-400 
wavenumber  region  with  decreasing  intensity  as  the  composition  changes  within  each  sample 
(figure  28). 

6.  MD  simulations  of  nanoindentation  (indenter  radius  =  9  nm)  into  fused  silica  were  conducted 
using  LAMMPS  (85)  using  a  pairwise  interatomic  potential  developed  using  novel  FM 
techniques  described  in  Izvekov  and  Rice  (14)  and  Izvekov  et  al.  (86).  However,  computed 
hardness  values  of  7.38  GPa  underestimated  experimental  values  determined  from  our  own 
indentation  experiments,  using  albeit  larger  indenters  (indenter  radius  =  3  fim)  (see  figure  30). 

7.  The  pairwise  interatomic  potential  that  was  developed  using  a  FM  technique  was  also  used  to 
determine  a  shock  Hugoniot  for  fused  silica  to  60  GPa;  it  was  initially  shown  to  be  in 
excellent  agreement  with  experimental  values  (see  figure  19  on  page  26  or  our  prior  year’s 
memorandum  report  (44).  Closer  examination  of  the  structure  of  our  silica  computational 
model  indicated  that  it  contained  regions  of  microcrystallinity  that  apparently  had  contributed 
to  our  excellent  prediction  of  the  shock  Hugoniot  for  fused  silica  (44).  Recomputation  of  the 
shock  Hugoniot  with  a  model  free  of  microcrystallites  resulted  in  a  Hugoniot  which 
overestimated  the  silica  density  relative  to  experimental  values  in  the  0-20  GPa  pressure 
range  (figure  35).  The  extension  of  the  pairwise  FM  model  to  include  angle  bending  forces, 
which  are  three-body  forces,  is  underway  to  improve  the  performance  of  the  model  at  low 
pressures  where  silica  densification  is  driven  by  structural  and  topological  changes.  At  this 
time,  however,  we  cannot  rule  out  entirely  the  importance  of  micro/nanocrystallinity  in  our 
modeling  of  glass,  since  work  by  Saito  et  al.  (32)  using  light  scattering  studies  concluded  that 
IRO  structures,  i.e.,  micro/nanocrystallites,  exist  in  silica  glass  (see  also  our  discussion  in 
section  4.2.7). 
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8.  First-principles  quantum  mechanical  methods  using  VASP  (62)  are  used  to  model 
densification  and  bulk  modulus  variations  with  pressure  in  fused  silica.  Ring  size 
distributions  (figure  46)  and  angle  distributions  for  Si  —  O  —  Si  and  O  —  Si  —  O  (figure  47) 
are  calculated,  and  future  work  will  consider  the  pressure  variation  of  these  distributions  and 
whether  they  can  be  validated  with  experimental  measurements  conducted  on  fused  silica 
under  extreme  pressure  using  diamond  anvil  cells.  FM  potentials  were  also  determined  for  a 
borosilicate  glass  system  (figure  52)  in  preparation  for  the  fiscal  year  2013  investigation  of  the 
role  of  boron  concentration  on  the  properties  of  glass  using  our  new  glass  processing  facility 
(see  section  13). 

9.  The  CASTEP  code  (63)  was  used  to  compute  initial  IR  and  Raman  spectra  for  both 
crystalline  and  amorphous  Si02.  It  was  found  that  the  calculated  IR  spectrum  (figure  53) 
reproduce  the  main  peaks  observed  in  experiments  on  a -quartz  (figure  54).  In  addition, 
calculated  Raman  spectra  (figure  60)  reproduced  the  three  main  regions  observed  in 
experiments  on  amorphous  silica  (figure  61). 

10.  Since  inelastic  deformation  and  a  polyamorphic  transformation  apparently  occurs 
simultaneously  in  some  glasses,  a  model  was  developed  to  account  for  this  coupled  behavior, 
which  captures  the  kinetics  of  the  polyamorphic  transformation  that  occur  in  plate  impact 
experiments  on  borosilicate  (compare  for  example  the  plate  impact  experiments  of  Alexander 
et  al.  (23)  in  figure  62a  with  simulation  results  in  figure  63b. 

1 1 .  We  have  conducted  a  dimensional  analysis  for  stable  crack  growth  in  brittle  elastic  solids 
subjected  to  both  indentation  and  wedge  loadings  and  found  that  for  Hertzian  cone  cracks 

D  =  (P/K )2/3,  where  D  is  the  width  of  the  base  of  the  cone  crack,  P  is  the  indentation  load, 
and  K  is  the  cohesive  modulus;  this  relation  is  modified  for  2-D  stress  and  deformation  fields 
l  =  (P/K)2,  where  l  is  the  crack  length.  This  last  similarity  relation  was  used  to  successfully 
validate  indentation  and  wedge-induced  crack  problems  using  the  peridynamics 
computational  code  that  was  developed  under  this  program.  Fully  3-D  deformation  and  stress 
fields  for  a  rigid  cylindrical  indenter  on  an  elastic  halfspace  were  presented,  which  are  useful 
for  validating  similar  elastic  field  predictions  by  computational  means  prior  to  densification 
and  failure  of  the  medium. 

12.  A  2-D  peridynamic  model  was  developed  to  model  nanoindentation  and  dynamic  crack 
propagation  in  fused  silica  as  an  extension  of  the  recent  1-D  model  of  Wildman  and 
Gazonas  (78);  a  PML  was  added  to  permit  infinite  computational  domains  that  are 
encountered  in  solution  to  certain  boundary  value  problems  (87).  Peridynamic  simulations  of 
2-D  cracking  induced  by  indentation  and  wedge  loadings  were  validated  using  the  universal 
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similarity  relation  /  =  ( PjK )2  derived  in  section  6.  A  3-D  axisymmetric  peridynamics 
computational  code  was  developed  and  verified  against  a  closed-form  solution  (25)  for 
vertical  displacements  in  an  elastic  halfspace  induced  by  a  flat,  cylindrical  indenter  illustrated 
in  figure  74. 

13.  Preliminary  material  characteristic  and  property  metrics:  At  this  point,  it  is  worthwhile  to  try 
articulating  a  set  of  material  characteristics  and  properties  that  may  have  significant  influence 
on  performance  of  silicate -based  glasses  in  various  ballistic  events  and  that  may  be  controlled 
by  systematic  chemical  substitutions.  Among  these  are  the  following:  linear  thermal 
expansion,  flash  thermal  conductivity,  inelastic  deformation  (plasticity),  actual  density  and 
calculated  free  volume,  macro-defects  (inclusions)  and  pores  (bubbles),  percentage  of 
micro/nanocry stalinity  and  nanochemical  inhomogeneities,  elastic  recovery,  and  controlled 
Hertzian  indentation.  The  quasi-static  and  dynamic  stress-strain  measurements  are  so 
susceptible  to  almost  unavoidable  surface  flaws  that  there  use  is  still  problematic. 

14.  The  current  glass  DSI  program  has  successfully  transitioned  into  a  core  mission  program 
within  WMRD  for  fiscal  year  2013  and  beyond.  An  important  focus  of  the  mission  program 
concerns  in-house  production  of  small  samples  of  glass  to  conduct  controlled 
composition- structure-property  relationship  studies  (see  section  13).  Chemical  analysis, 
quasi-static  property  characterization,  dynamic  properties  characterization,  and  ballistic 
performance  will  be  measured  and  reported  on  a  variety  of  chemically  substituted  glasses; 
data  and  measurements  from  this  effort  will  provide  a  means  to  validate  a  parallel 
computational  effort  using  DFT  and  MD  methods  for  the  “design”  of  ballistically  enhanced 
glass  formulations,  which  have  not  yet  been  synthesized! 
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non-crystalline  solids,  International  Congress  on  Glass,  Prague,  1-5  July  2013 

'Paper  accepted  yet  withdrawn  for  apparent  budgetary  considerations 
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15.2  Publications 


Wildman,  R.  A.;  Gazonas,  G.  A.  A  perfectly  matched  layer  for  peridynamics  in  one 
dimension ;  ARL-TR-5626;  U.S.  Army  Research  Laboratory:  Aberdeen  Proving  Ground, 
MD,  August,  2011. 

Wildman,  R.  A.;  Gazonas,  G.  A.  A  perfectly  matched  layer  for  peridynamics  in  two 
dimensions.  J.Mech.  Mater,  and  Struct.  2013,  7  (8-9),  765-781. 

Izvekov,  S.;  Rice,  B.  M.  Mechanism  of  densification  in  silica  glass  under  pressure  as  revealed 
by  bottom-up  pairwise  effective  interaction  model;  J.  Chem.  Phys.  2012,  136  (13),  134508. 

*Grujicic,  M.;  Bell,  W.  C.;  Pandurangan,  B.;  Cheeseman,  B.A.;  Patel,  P.;  Gazonas,  G.A. 
Inclusion  of  material  non-linearity  and  inelasticity  into  a  continuum-level  material  model  for 
soda-lime  glass.  J.  Mater,  and Des.  2012,35,144-155. 

*  Souza,  F.;  Allen,  D.  H.;  Modeling  the  transition  of  microcracks  into  macrocracks  in 
heterogeneous  media  using  a  two-way  coupled  multiscale  model.  Inti.  J.  Solids  Structures 
2011,  48,  3160-3175.  (Scientific  Services  Program  administered  by  Battelle  under  contract 
W9 1 1 NF-07  -D  -000 1 ;  TCN  09-201). 

*Grujicic,  M.;  Pandurangan,  B.;  Zhang,  Z.;  Bell,  W.C.;  Gazonas,  G.  A.;  Patel,  P;  Cheeseman, 
B.  A.  Molecular-level  analysis  of  shock-wave  physics  and  derivation  of  the  Hugoniot  relations 
for  fused  silica.  J.  Mater.  Eng.  2011,  Perf.,  DOI:  10. 1007/s  1 1665-01 1-0005-2. 

*Grujicic,  M.;  Pandurangan,  B.;  Bell,  W.  C.;  Cheeseman,  B.  A.;  Patel,  P;  Gazonas,  G.  A. 
Molecular-level  analysis  of  shock- wave  physics  and  derivation  of  the  Hugoniot  relations  for 
soda-lime  glass.  J.  of  Mater.  Sci.  2011,  46  (22),  7298-7312. 

Gazonas,  G.  A.;  McCauley,  J.  W.;  Batyrev,  I.  G.;  Becker,  R.  C.;  Izvekov,  S.;  Jenkins,  T.  A.; 
Patel,  P;  Rice,  B.  M.;  Schuster,  B.  E.;  Weingarten,  N.  S.;  Wildman,  R.  A.  Multiscale  modeling 
of  non-crystalline  ceramics  (glass)  (final  report)-,  ARL-TR-6353;  U.S.  Army  Research 
Laboratory:  Aberdeen  Proving  Ground,  MD,  February  2013. 

Gazonas,  G.  A.;  McCauley,  J.  W.;  Batyrev,  I.  G.;  Becker,  R.  C.;  Izvekov,  S.;  Patel,  P;  Rice,  B. 
M.;  Schuster,  B.  E.;  Weingarten,  N.  S.;  Wildman,  R.  A.  Multiscale  modeling  of 
non-crystalline  ceramics  (glass)  ( FY11 );  ARL-MR-0802;  U.S.  Army  Research  Laboratory: 
Aberdeen  Proving  Ground,  MD,  January  2012. 

^Funded  in  part  by  the  Glass  DSI  program. 
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Gazonas,  G.  A.;  McCauley,  J.  W.;  Batyrev,  I.  G.;  Becker,  R.  C.;  Patel,  R;  Rice,  B.  M.; 
Weingarten,  N.S.  Multiscale  modeling  of  non- crystalline  ceramics  (glass)',  ARL-MR-0765; 
U.S.  Army  Research  Laboratory:  Aberdeen  Proving  Ground,  MD,  February  201 1 . 

Strassburger,  E.;  Bauer,  S.;  Hunzinger,  M.;  Analysis  of  the  fragmentation  ofAlONand 
bi-modal  grain  sized  Spinel,  Fined  Report  Report  E  36/12',  Research  period: 
10/2010-10/2011;  Contract  No.  W91  INF-10-2-0100  Contractor:  Fraunhofer-Gesellschaft 
zur  Forderung  der  angewandten  Forschung  e.V. 

■^Becker,  R.  Formulation  of  a  glass  model  to  capture  observations  from  high-rate  ballistic 
penetration',  ARL-TR-6086;  U.S.  Army  Research  Laboratory:  Aberdeen  Proving  Ground, 
MD,  August,  2012. 

15.3  Personnel  Hires 

An  Oak  Ridge  Institute  of  Science  and  Education  (ORISE)  post-doctoral  fellow  was  hired  in 
FY12  and  offered  a  full-time  government  position  in  FY13.  An  ORISE  summer  student  was 
hired  in  FY 1 1 . 

15.4  Transition  to  Core 

The  DSI  program  has  been  transitioned  to  a  core  mission  program  in  WMRD  for  fiscal  year 
2013  and  beyond.  In  addition,  a  glass  processing  facility  (infrastructure  enhancement)  has 
been  established  in  WMRD  to  produce  small  samples  of  glass  to  conduct 
composition-structure-property  relationship  studies  as  part  of  the  new  core  mission  program. 

15.5  Computational  Code  Development 

The  1-,  2-,  and  3-D  axisymmetric  peridynamics  codes  with  PMLs  capable  of  modeling  wave 
propagation  and  fracture  in  brittle  materials  have  been  developed  and  verified  under  this  DSI 
program. 

15.6  Other 

A  one-day  short  course  on  the  “Fundamentals  of  Glass  Science,”  which  was  taught  by 
Professor  Arun  K.  Varshneya,  Alfred  University,  at  ARL  on  October  29,  2010.  Course 
attendees  received  copies  of  Varshneya’s  Fundamentals  of  Inorganic  Glasses  (26). 


96 


A  short-term  conceptual  project  to  determine  an  effective  experimental  and  theoretical 
approach  to  model  and  characterize  the  role  of  glassy  materials  in  resisting  ballistic  impact 
was  conducted  by  Professor  Richard  Lehman,  Rutgers  University.  A  short  synopsis  of  this 
project  was  briefed  at  ARL  on  17  November,  201 1  (section  12). 

Seminars:  “Atomistic  Simulations  of  Vitreous  Silica,”  Cormack,  A.  N.,  and  Inamori,  K., 
School  of  Engineering,  NY  State  College  of  Ceramics,  Alfred  University,  9  January  2012. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


1-D 

one-dimensional 

2-D 

two-dimensional 

a-SiO-2 

crystalline  quartz 

Al 

aluminum 

AlON 

aluminum  oxynitride 

ARL 

U.S.  Army  Research  Laboratory 

AIMD 

ab  initio  molecular  dynamics 

a-SiC>2 

fused  silica  or  amorphous  quartz 

B 

boron 

Ca 

calcium 

CSM 

continuous  stiffness  measurement 

DAC 

diamond  anvil  cell 

DSI 

Director’s  Strategic  Initiative 

EOI 

edge-on-impact 

EOS 

equation  of  state 

FEM 

finite  element  method 

FIB 

focused-ion-beam 

HEL 

Hugoniot  elastic  limit 

IRO 

intermediate -range  order 

LAMMPS 

large-scale  atomic/molecular  massively  parallel  simulator 

LB 

less  brittle 

MD 

molecular  dynamics 

Na 

sodium 

MEDE 

materials  in  extreme  dynamic  environments 

0 

oxygen 

PAW 

projector  augmented  wave 

PML 

perfectly  matched  layer 

RDF 

radial  distribution  function 

SAXS 

small-angle  x-ray  scattering 

SCJ 

shaped-charge  jet 

SEET 

strain  energy  equivalence  theory 

SEM 

scanning  electron  microscopy 
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Si 

silicon 

SiC 

silicon  carbide 

SL 

soda  lime 

SRO 

short-range  order 

STEM 

scanning  tunneling  electron  microscopy 

VASP 

Vienna  ab  initio  simulation  package 

WMRD 

Weapons  and  Materials  Research  Directorate 
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